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 // AliTPCCorrection class //
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 //
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. //
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 //
33 // The class allows "effective Omega Tau" corrections to be shifted to the //
34 // single distortion classes. //
36 // Note: This class is normally used via the class AliTPCComposedCorrection //
38 // date: 27/04/2010 //
39 // Authors: Magnus Mager, Stefan Rossegger, Jim Thomas //
40 ////////////////////////////////////////////////////////////////////////////////
41 #include "Riostream.h"
46 #include <TTreeStream.h>
49 #include <TTimeStamp.h>
50 #include <AliCDBStorage.h>
52 #include <AliCDBMetaData.h>
54 #include "AliExternalTrackParam.h"
55 #include "AliTrackPointArray.h"
56 #include "TDatabasePDG.h"
57 #include "AliTrackerBase.h"
58 #include "AliTPCROC.h"
61 #include "AliTPCCorrection.h"
63 ClassImp(AliTPCCorrection)
65 // FIXME: the following values should come from the database
66 const Double_t AliTPCCorrection::fgkTPC_Z0 =249.7; // nominal gating grid position
67 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
68 const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
69 const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
70 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
71 const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
74 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
75 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
76 84.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
77 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
78 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
79 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
80 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
81 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
82 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
83 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
84 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
85 249.0, 249.5, 250.0, 251.5, 252.0 } ;
87 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
88 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
89 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
90 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
91 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
92 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
93 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
94 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
95 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
96 -100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
97 -80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
98 -60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
99 -40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
100 -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
101 -1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
102 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
103 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
104 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
105 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
106 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
107 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
108 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
109 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
110 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
111 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
112 202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
113 222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
114 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
118 AliTPCCorrection::AliTPCCorrection()
119 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
122 // default constructor
126 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
127 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
130 // default constructor, that set the name and title
134 AliTPCCorrection::~AliTPCCorrection() {
136 // virtual destructor
140 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
142 // Corrects the initial coordinates x (cartesian coordinates)
143 // according to the given effect (inherited classes)
144 // roc represents the TPC read out chamber (offline numbering convention)
147 GetCorrection(x,roc,dx);
148 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
151 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
153 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
154 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
155 // roc represents the TPC read out chamber (offline numbering convention)
158 GetCorrection(x,roc,dx);
159 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
162 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
164 // Distorts the initial coordinates x (cartesian coordinates)
165 // according to the given effect (inherited classes)
166 // roc represents the TPC read out chamber (offline numbering convention)
169 GetDistortion(x,roc,dx);
170 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
173 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
175 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
176 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
177 // roc represents the TPC read out chamber (offline numbering convention)
180 GetDistortion(x,roc,dx);
181 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
184 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
186 // This function delivers the correction values dx in respect to the inital coordinates x
187 // roc represents the TPC read out chamber (offline numbering convention)
188 // Note: The dx is overwritten by the inherited effectice class ...
190 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
193 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
195 // This function delivers the distortion values dx in respect to the inital coordinates x
196 // roc represents the TPC read out chamber (offline numbering convention)
198 GetCorrection(x,roc,dx);
199 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
202 void AliTPCCorrection::Init() {
204 // Initialization funtion (not used at the moment)
208 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
214 void AliTPCCorrection::Print(Option_t* /*option*/) const {
216 // Print function to check which correction classes are used
217 // option=="d" prints details regarding the setted magnitude
218 // option=="a" prints the C0 and C1 coefficents for calibration purposes
220 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
223 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
225 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
226 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
231 //SetOmegaTauT1T2(omegaTau, t1, t2);
234 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
236 // Simple plot functionality.
237 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
238 // in respect to position z within the XY plane.
239 // The histogramm has nx times ny entries.
242 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
243 nx,-250.,250.,ny,-250.,250.);
246 Int_t roc=z>0.?0:18; // FIXME
247 for (Int_t iy=1;iy<=ny;++iy) {
248 x[1]=h->GetYaxis()->GetBinCenter(iy);
249 for (Int_t ix=1;ix<=nx;++ix) {
250 x[0]=h->GetXaxis()->GetBinCenter(ix);
251 GetCorrection(x,roc,dx);
252 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
253 if (90.<=r0 && r0<=250.) {
254 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
255 h->SetBinContent(ix,iy,r1-r0);
258 h->SetBinContent(ix,iy,0.);
264 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
266 // Simple plot functionality.
267 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
268 // in respect to position z within the XY plane.
269 // The histogramm has nx times ny entries.
272 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
273 nx,-250.,250.,ny,-250.,250.);
276 Int_t roc=z>0.?0:18; // FIXME
277 for (Int_t iy=1;iy<=ny;++iy) {
278 x[1]=h->GetYaxis()->GetBinCenter(iy);
279 for (Int_t ix=1;ix<=nx;++ix) {
280 x[0]=h->GetXaxis()->GetBinCenter(ix);
281 GetCorrection(x,roc,dx);
282 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
283 if (90.<=r0 && r0<=250.) {
284 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
285 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
287 Float_t dphi=phi1-phi0;
288 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
289 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
291 h->SetBinContent(ix,iy,r0*dphi);
294 h->SetBinContent(ix,iy,0.);
300 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
302 // Simple plot functionality.
303 // Returns a 2d hisogram which represents the corrections in r direction (dr)
304 // in respect to angle phi within the ZR plane.
305 // The histogramm has nx times ny entries.
307 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
308 nz,-250.,250.,nr,85.,250.);
310 for (Int_t ir=1;ir<=nr;++ir) {
311 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
312 x[0]=radius*TMath::Cos(phi);
313 x[1]=radius*TMath::Sin(phi);
314 for (Int_t iz=1;iz<=nz;++iz) {
315 x[2]=h->GetXaxis()->GetBinCenter(iz);
316 Int_t roc=x[2]>0.?0:18; // FIXME
317 GetCorrection(x,roc,dx);
318 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
319 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
320 h->SetBinContent(iz,ir,r1-r0);
328 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
330 // Simple plot functionality.
331 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
332 // in respect to angle phi within the ZR plane.
333 // The histogramm has nx times ny entries.
335 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
336 nz,-250.,250.,nr,85.,250.);
338 for (Int_t iz=1;iz<=nz;++iz) {
339 x[2]=h->GetXaxis()->GetBinCenter(iz);
340 Int_t roc=x[2]>0.?0:18; // FIXME
341 for (Int_t ir=1;ir<=nr;++ir) {
342 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
343 x[0]=radius*TMath::Cos(phi);
344 x[1]=radius*TMath::Sin(phi);
345 GetCorrection(x,roc,dx);
346 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
347 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
348 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
350 Float_t dphi=phi1-phi0;
351 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
352 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
354 h->SetBinContent(iz,ir,r0*dphi);
360 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
361 const char *xlabel,const char *ylabel,const char *zlabel,
362 Int_t nbinsx,Double_t xlow,Double_t xup,
363 Int_t nbinsy,Double_t ylow,Double_t yup) {
365 // Helper function to create a 2d histogramm of given size
371 while (gDirectory->FindObject(hname.Data())) {
378 TH2F *h=new TH2F(hname.Data(),title,
381 h->GetXaxis()->SetTitle(xlabel);
382 h->GetYaxis()->SetTitle(ylabel);
383 h->GetZaxis()->SetTitle(zlabel);
389 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
391 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
392 const Double_t er[kNZ][kNR], Double_t &er_value )
395 // Interpolate table - 2D interpolation
397 Double_t save_er[10] ;
399 Search( kNZ, fgkZList, z, fJLow ) ;
400 Search( kNR, fgkRList, r, fKLow ) ;
401 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
402 if ( fKLow < 0 ) fKLow = 0 ;
403 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
404 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
406 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
407 save_er[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
409 er_value = Interpolate( &fgkZList[fJLow], save_er, order, z ) ;
414 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
415 const Int_t order, const Double_t x )
418 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
422 if ( order == 2 ) { // Quadratic Interpolation = 2
423 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
424 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
425 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
426 } else { // Linear Interpolation = 1
427 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
435 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
438 // Search an ordered table by starting at the most recently used point
441 Long_t middle, high ;
442 Int_t ascend = 0, increment = 1 ;
444 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
446 if ( low < 0 || low > n-1 ) {
447 low = -1 ; high = n ;
448 } else { // Ordered Search phase
449 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
450 if ( low == n-1 ) return ;
452 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
455 high = low + increment ;
456 if ( high > n-1 ) { high = n ; break ; }
459 if ( low == 0 ) { low = -1 ; return ; }
461 while ( (Int_t)( x < xArray[low] ) == ascend ) {
464 if ( increment >= high ) { low = -1 ; break ; }
465 else low = high - increment ;
470 while ( (high-low) != 1 ) { // Binary Search Phase
471 middle = ( high + low ) / 2 ;
472 if ( (Int_t)( x >= xArray[middle] ) == ascend )
478 if ( x == xArray[n-1] ) low = n-2 ;
479 if ( x == xArray[0] ) low = 0 ;
484 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){
486 // Fit the track parameters - without and with distortion
487 // 1. Space points in the TPC are simulated along the trajectory
488 // 2. Space points distorted
489 // 3. Fits the non distorted and distroted track to the reference plane at refX
490 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
492 // trackIn - input track parameters
493 // refX - reference X to fit the track
494 // dir - direction - out=1 or in=-1
495 // pcstream - debug streamer to check the results
497 // see AliExternalTrackParam.h documentation:
498 // track1.fP[0] - local y (rphi)
500 // track1.fP[2] - sinus of local inclination angle
501 // track1.fP[3] - tangent of deep angle
502 // track1.fP[4] - 1/pt
503 AliTPCROC * roc = AliTPCROC::Instance();
504 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
505 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
506 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
508 const Double_t kMaxSnp = 0.85;
509 const Double_t kSigmaY=0.1;
510 const Double_t kSigmaZ=0.1;
511 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
513 AliExternalTrackParam track(trackIn); //
515 AliTrackPointArray pointArray0(npoints0);
516 AliTrackPointArray pointArray1(npoints0);
518 AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
520 // simulate the track
522 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
523 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
524 AliTrackerBase::PropagateTrackTo(&track,radius,kMass,3,kTRUE,kMaxSnp);
526 AliTrackPoint pIn0; // space point
528 Int_t sector= (xyz[2]>0)? 0:18;
529 pointArray0.GetPoint(pIn0,npoints);
530 pointArray1.GetPoint(pIn1,npoints);
531 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
532 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
533 DistortPoint(distPoint, sector);
534 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
535 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
538 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
539 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
540 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
541 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
542 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
543 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
544 pointArray0.AddPoint(npoints, &pIn0);
545 pointArray1.AddPoint(npoints, &pIn1);
547 if (npoints>=npoints0) break;
552 AliExternalTrackParam *track0=0;
553 AliExternalTrackParam *track1=0;
554 AliTrackPoint point1,point2,point3;
555 if (dir==1) { //make seed inner
556 pointArray0.GetPoint(point1,1);
557 pointArray0.GetPoint(point2,10);
558 pointArray0.GetPoint(point3,20);
560 if (dir==-1){ //make seed outer
561 pointArray0.GetPoint(point1,npoints-20);
562 pointArray0.GetPoint(point2,npoints-10);
563 pointArray0.GetPoint(point3,npoints-1);
565 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
566 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
569 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
570 Int_t ipoint= (dir>0) ? ipoint: npoints-1-jpoint;
574 pointArray0.GetPoint(pIn0,ipoint);
575 pointArray1.GetPoint(pIn1,ipoint);
576 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
577 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
579 AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,1,kFALSE,kMaxSnp);
580 AliTrackerBase::PropagateTrackTo(track1,prot1.GetX(),kMass,1,kFALSE,kMaxSnp);
583 Double_t pointPos[2]={0,0};
584 Double_t pointCov[3]={0,0,0};
585 pointPos[0]=prot0.GetY();//local y
586 pointPos[1]=prot0.GetZ();//local z
587 pointCov[0]=prot0.GetCov()[3];//simay^2
588 pointCov[1]=prot0.GetCov()[4];//sigmayz
589 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
590 track0->Update(pointPos,pointCov);
592 pointPos[0]=prot1.GetY();//local y
593 pointPos[1]=prot1.GetZ();//local z
594 pointCov[0]=prot1.GetCov()[3];//simay^2
595 pointCov[1]=prot1.GetCov()[4];//sigmayz
596 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
597 track1->Update(pointPos,pointCov);
600 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,2.,kTRUE,kMaxSnp);
601 track1->Rotate(track0->GetAlpha());
602 track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
604 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
605 "point0.="<<&pointArray0<< // points
606 "point1.="<<&pointArray1<< // distorted points
607 "trackIn.="<<&track<< // original track
608 "track0.="<<track0<< // fitted track
609 "track1.="<<track1<< // fitted distorted track
611 new(&trackIn) AliExternalTrackParam(*track0);
620 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
622 // create the distortion tree on a mesh with granularity given by step
623 // return the tree with distortions at given position
624 // Map is created on the mesh with given step size
626 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
628 for (Double_t x= -250; x<250; x+=step){
629 for (Double_t y= -250; y<250; y+=step){
630 Double_t r = TMath::Sqrt(x*x+y*y);
633 for (Double_t z= -250; z<250; z+=step){
634 Int_t roc=(z>0)?0:18;
638 Double_t phi = TMath::ATan2(y,x);
639 DistortPoint(xyz,roc);
640 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
641 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
642 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
643 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
644 Double_t dx = xyz[0]-x;
645 Double_t dy = xyz[1]-y;
646 Double_t dz = xyz[2]-z;
648 Double_t drphi=(phi1-phi)*r;
649 (*pcstream)<<"distortion"<<
650 "x="<<x<< // original position
655 "x1="<<xyz[0]<< // distorted position
661 "dx="<<dx<< // delta position
671 TFile f(Form("correction%s.root",GetName()));
672 TTree * tree = (TTree*)f.Get("distortion");
673 TTree * tree2= tree->CopyTree("1");
674 tree2->SetName(Form("dist%s",GetName()));
675 tree2->SetDirectory(0);
683 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, TObjArray * corrArray, Int_t step, Bool_t debug ){
686 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
687 // calculates partial distortions
688 // Partial distortion is stored in the resulting tree
689 // Output is storred in the file distortion_<dettype>_<partype>.root
690 // Partial distortion is stored with the name given by correction name
693 // Parameters of function:
694 // input - input tree
695 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
696 // ppype - parameter type
697 // corrArray - array with partial corrections
698 // step - skipe entries - if 1 all entries processed - it is slow
699 // debug 0 if debug on also space points dumped - it is slow
700 const Int_t kMinEntries=50;
701 Double_t phi,theta, snp, mean,rms, entries;
702 tinput->SetBranchAddress("theta",&theta);
703 tinput->SetBranchAddress("phi", &phi);
704 tinput->SetBranchAddress("snp",&snp);
705 tinput->SetBranchAddress("mean",&mean);
706 tinput->SetBranchAddress("rms",&rms);
707 tinput->SetBranchAddress("entries",&entries);
708 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
710 Int_t nentries=tinput->GetEntries();
711 Int_t ncorr=corrArray->GetEntries();
712 Double_t corrections[100]; //
714 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
717 if (dtype==0) {refX=85; dir=-1;}
718 if (dtype==1) {refX=245; dir=1;}
719 if (dtype==2) {refX=0; dir=-1;}
721 for (Int_t ientry=0; ientry<nentries; ientry+=step){
722 tinput->GetEntry(ientry);
727 tPar[4]=0.00001; // should be calculated - non equal to 0
729 cout<<"Entry\n\n"<<ientry<<endl;
730 cout<<"dtype="<<dtype<< // detector match type
731 "ptype="<<ptype<< // parameter type
732 "theta="<<theta<< // theta
735 "mean="<<mean<< // mean dist value
737 "entries="<<entries<<endl; // number of entries in bin
739 if (TMath::Abs(snp)>0.251) continue;
741 "dtype="<<dtype<< // detector match type
742 "ptype="<<ptype<< // parameter type
743 "theta="<<theta<< // theta
746 "mean="<<mean<< // mean dist value
748 "entries="<<entries;// number of entries in bin
750 for (Int_t icorr=0; icorr<ncorr; icorr++) {
751 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
752 corrections[icorr]=0;
753 if (entries>kMinEntries){
754 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
755 AliExternalTrackParam *trackOut = 0;
756 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
757 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
758 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
762 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
764 (*pcstream)<<"fit"<<"\n";
772 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
774 // Store object in the OCDB
775 // By default the object is stored in the current directory
776 // default comment consit of user name and the date
778 TString ocdbStorage="";
779 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
780 AliCDBMetaData *metaData= new AliCDBMetaData();
781 metaData->SetObjectClassName("AliTPCCorrection");
782 metaData->SetResponsible("Marian Ivanov");
783 metaData->SetBeamPeriod(1);
784 metaData->SetAliRootVersion("05-25-01"); //root version
785 TString userName=gSystem->GetFromPipe("echo $USER");
786 TString date=gSystem->GetFromPipe("date");
788 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
789 if (comment) metaData->SetComment(comment);
791 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
792 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
793 gStorage->Put(this, (*id1), metaData);