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>
56 #include "AliTPCTransform.h"
57 #include "AliTPCcalibDB.h"
58 #include "AliTPCExB.h"
59 #include "AliTPCCorrection.h"
60 #include "AliTPCRecoParam.h"
62 #include "AliExternalTrackParam.h"
63 #include "AliTrackPointArray.h"
64 #include "TDatabasePDG.h"
65 #include "AliTrackerBase.h"
66 #include "AliTPCROC.h"
67 #include "THnSparse.h"
68 #include "AliTPCLaserTrack.h"
70 #include "AliTPCCorrection.h"
72 ClassImp(AliTPCCorrection)
74 // FIXME: the following values should come from the database
75 const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
76 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
77 const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
78 const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
79 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
80 const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
83 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
84 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
85 84.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
86 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
87 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
88 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
89 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
90 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
91 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
92 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
93 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
94 249.0, 249.5, 250.0, 251.5, 252.0 } ;
96 const 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,
116 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
117 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
118 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
119 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
120 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
121 202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
122 222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
123 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
127 AliTPCCorrection::AliTPCCorrection()
128 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
131 // default constructor
135 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
136 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
139 // default constructor, that set the name and title
143 AliTPCCorrection::~AliTPCCorrection() {
145 // virtual destructor
149 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
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)
156 GetCorrection(x,roc,dx);
157 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
160 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
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)
167 GetCorrection(x,roc,dx);
168 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
171 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
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)
178 GetDistortion(x,roc,dx);
179 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
182 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
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)
189 GetDistortion(x,roc,dx);
190 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
193 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
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 ...
199 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
202 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
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)
207 GetCorrection(x,roc,dx);
208 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
211 void AliTPCCorrection::Init() {
213 // Initialization funtion (not used at the moment)
217 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
223 void AliTPCCorrection::Print(Option_t* /*option*/) const {
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
229 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
232 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
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
240 //SetOmegaTauT1T2(omegaTau, t1, t2);
243 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
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.
251 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
252 nx,-250.,250.,ny,-250.,250.);
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);
267 h->SetBinContent(ix,iy,0.);
273 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
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.
281 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
282 nx,-250.,250.,ny,-250.,250.);
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]);
296 Float_t dphi=phi1-phi0;
297 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
298 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
300 h->SetBinContent(ix,iy,r0*dphi);
303 h->SetBinContent(ix,iy,0.);
309 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
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.
316 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
317 nz,-250.,250.,nr,85.,250.);
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);
337 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
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.
344 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
345 nz,-250.,250.,nr,85.,250.);
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]);
359 Float_t dphi=phi1-phi0;
360 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
361 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
363 h->SetBinContent(iz,ir,r0*dphi);
369 TH2F* 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) {
374 // Helper function to create a 2d histogramm of given size
380 while (gDirectory->FindObject(hname.Data())) {
387 TH2F *h=new TH2F(hname.Data(),title,
390 h->GetXaxis()->SetTitle(xlabel);
391 h->GetYaxis()->SetTitle(ylabel);
392 h->GetZaxis()->SetTitle(zlabel);
398 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
400 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
401 const Double_t er[kNZ][kNR], Double_t &erValue ) {
403 // Interpolate table - 2D interpolation
405 Double_t saveEr[10] ;
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 ;
414 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
415 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
417 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
422 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
423 const Int_t order, const Double_t x ) {
425 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
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] ) ;
442 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
444 // Search an ordered table by starting at the most recently used point
447 Long_t middle, high ;
448 Int_t ascend = 0, increment = 1 ;
450 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
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 ;
458 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
461 high = low + increment ;
462 if ( high > n-1 ) { high = n ; break ; }
465 if ( low == 0 ) { low = -1 ; return ; }
467 while ( (Int_t)( x < xArray[low] ) == ascend ) {
470 if ( increment >= high ) { low = -1 ; break ; }
471 else low = high - increment ;
476 while ( (high-low) != 1 ) { // Binary Search Phase
477 middle = ( high + low ) / 2 ;
478 if ( (Int_t)( x >= xArray[middle] ) == ascend )
484 if ( x == xArray[n-1] ) low = n-2 ;
485 if ( x == xArray[0] ) low = 0 ;
490 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
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
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
503 // see AliExternalTrackParam.h documentation:
504 // track1.fP[0] - local y (rphi)
506 // track1.fP[2] - sinus of local inclination angle
507 // track1.fP[3] - tangent of deep angle
508 // track1.fP[4] - 1/pt
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);
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();
519 AliExternalTrackParam track(trackIn); //
521 AliTrackPointArray pointArray0(npoints0);
522 AliTrackPointArray pointArray1(npoints0);
524 AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
526 // simulate the track
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++){
530 AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
532 xyz[0]+=gRandom->Gaus(0,0.005);
533 xyz[1]+=gRandom->Gaus(0,0.005);
534 xyz[2]+=gRandom->Gaus(0,0.005);
535 AliTrackPoint pIn0; // space point
537 Int_t sector= (xyz[2]>0)? 0:18;
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]};
542 DistortPoint(distPoint, sector);
543 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
544 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
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);
556 if (npoints>=npoints0) break;
558 if (npoints<npoints0/2) return 0;
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);
567 pointArray0.GetPoint(point2,30);
568 pointArray0.GetPoint(point3,60);
570 if (dir==-1){ //make seed outer
571 pointArray0.GetPoint(point1,npoints-60);
572 pointArray0.GetPoint(point2,npoints-30);
573 pointArray0.GetPoint(point3,npoints-1);
575 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
576 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
578 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
579 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
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
588 AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
589 AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
590 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
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);
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
605 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
606 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
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);
613 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
614 track1->Rotate(track0->GetAlpha());
615 track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
617 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
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
624 new(&trackIn) AliExternalTrackParam(*track0);
633 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
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
639 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
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);
646 for (Double_t z= -250; z<250; z+=step){
647 Int_t roc=(z>0)?0:18;
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;
661 Double_t drphi=(phi1-phi)*r;
662 (*pcstream)<<"distortion"<<
663 "x="<<x<< // original position
668 "x1="<<xyz[0]<< // distorted position
674 "dx="<<dx<< // delta position
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);
696 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
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
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
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;
716 const Int_t kMinEntries=50;
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));
726 Int_t nentries=tinput->GetEntries();
727 Int_t ncorr=corrArray->GetEntries();
728 Double_t corrections[100]={0}; //
730 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
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;}
738 for (Int_t ientry=0; ientry<nentries; ientry+=step){
739 tinput->GetEntry(ientry);
740 if (TMath::Abs(snp)>kMaxSnp) continue;
745 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
746 Double_t bz=AliTrackerBase::GetBz();
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;
749 AliExternalTrackParam track(refX,phi,tPar,cov);
753 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
754 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
756 "bz="<<bz<< // magnetic filed used
757 "dtype="<<dtype<< // detector match type
758 "ptype="<<ptype<< // parameter type
759 "theta="<<theta<< // theta
762 "mean="<<mean<< // mean dist value
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
769 "entries="<<entries;// number of entries in bin
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);
779 if (dtype==0) {refX=85.; dir=-1;}
780 if (dtype==1) {refX=275.; dir=1;}
781 if (dtype==2) {refX=0; dir=-1;}
782 if (dtype==3) {refX=360.; dir=-1;}
785 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
786 trackOut->Rotate(trackIn.GetAlpha());
787 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
789 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
792 corrections[icorr]=0;
794 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
798 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
799 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
800 // for points it is neccessary
802 (*pcstream)<<"fit"<<"\n";
809 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
811 // Make a laser fit tree for global minimization
813 const Double_t cutErrY=0.1;
814 const Double_t cutErrZ=0.1;
815 const Double_t kEpsilon=0.00000001;
820 AliTPCLaserTrack *ltr=0;
821 AliTPCLaserTrack::LoadTracks();
822 tree->SetBranchAddress("dY.",&vecdY);
823 tree->SetBranchAddress("dZ.",&vecdZ);
824 tree->SetBranchAddress("eY.",&veceY);
825 tree->SetBranchAddress("eZ.",&veceZ);
826 tree->SetBranchAddress("LTr.",<r);
827 Int_t entries= tree->GetEntries();
828 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
829 Double_t bz=AliTrackerBase::GetBz();
832 for (Int_t ientry=0; ientry<entries; ientry++){
833 tree->GetEntry(ientry);
834 if (!ltr->GetVecGX()){
837 TVectorD * delta= (itype==0)? vecdY:vecdZ;
838 TVectorD * err= (itype==0)? veceY:veceZ;
840 for (Int_t irow=0; irow<159; irow++){
841 Int_t nentries = 1000;
842 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
843 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
845 Double_t phi =(*ltr->GetVecPhi())[irow];
846 Double_t theta =ltr->GetTgl();
847 Double_t mean=delta->GetMatrixArray()[irow];
848 Double_t gx=0,gy=0,gz=0;
849 Double_t snp = (*ltr->GetVecP2())[irow];
850 Double_t rms = 0.1+err->GetMatrixArray()[irow];
851 gx = (*ltr->GetVecGX())[irow];
852 gy = (*ltr->GetVecGY())[irow];
853 gz = (*ltr->GetVecGZ())[irow];
854 Int_t bundle= ltr->GetBundle();
857 // get delta R used in reconstruction
858 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
859 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
860 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
861 Double_t xyz0[3]={gx,gy,gz};
862 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
864 // old ExB correction
866 if(recoParam&&recoParam->GetUseExBCorrection()) {
867 Double_t xyz1[3]={gx,gy,gz};
868 calib->GetExB()->Correct(xyz0,xyz1);
869 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
872 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
873 Float_t xyz1[3]={gx,gy,gz};
874 Int_t sector=(gz>0)?0:18;
875 correction->CorrectPoint(xyz1, sector);
876 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
882 "bz="<<bz<< // magnetic filed used
883 "dtype="<<dtype<< // detector match type
884 "ptype="<<itype<< // parameter type
885 "theta="<<theta<< // theta
888 "mean="<<mean<< // mean dist value
890 "gx="<<gx<< // global position
891 "gy="<<gy<< // global position
892 "gz="<<gz<< // global position
893 "dRrec="<<dRrec<< // delta Radius in reconstruction
894 "id="<<bundle<< //bundle
895 "entries="<<nentries;// number of entries in bin
898 Double_t ky = TMath::Tan(TMath::ASin(snp));
899 Int_t ncorr = corrArray->GetEntries();
900 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
901 Double_t phi0 = TMath::ATan2(gy,gx);
902 Double_t distortions[1000]={0};
903 Double_t distortionsR[1000]={0};
904 for (Int_t icorr=0; icorr<ncorr; icorr++) {
905 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
906 Float_t distPoint[3]={gx,gy,gz};
907 Int_t sector= (gz>0)? 0:18;
909 corr->DistortPoint(distPoint, sector);
911 Double_t value=distPoint[2]-gz;
913 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
914 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
915 Double_t drphi= r0*(phi1-phi0);
917 distortions[icorr] = drphi-ky*dr;
918 distortionsR[icorr] = dr;
921 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
922 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
924 (*pcstream)<<"fit"<<"\n";
932 void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
934 // make a distortion map out ou fthe residual histogram
935 // Results are written to the debug streamer - pcstream
937 // his0 - input (4D) residual histogram
938 // pcstream - file to write the tree
940 // marian.ivanov@cern.ch
941 const Int_t kMinEntries=50;
942 Int_t nbins1=his0->GetAxis(1)->GetNbins();
943 Int_t first1=his0->GetAxis(1)->GetFirst();
944 Int_t last1 =his0->GetAxis(1)->GetLast();
946 Double_t bz=AliTrackerBase::GetBz();
947 Int_t idim[4]={0,1,2,3};
948 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
950 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
951 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
952 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
953 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
954 Int_t first3 = his1->GetAxis(3)->GetFirst();
955 Int_t last3 = his1->GetAxis(3)->GetLast();
958 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
959 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
960 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
962 his1->GetAxis(3)->SetRangeUser(-1,1);
965 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
966 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
967 Int_t first2 = his3->GetAxis(2)->GetFirst();
968 Int_t last2 = his3->GetAxis(2)->GetLast();
970 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
971 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
972 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
973 TH1 * hisDelta = his3->Projection(0);
975 Double_t entries = hisDelta->GetEntries();
976 Double_t mean=0, rms=0;
977 if (entries>kMinEntries){
978 mean = hisDelta->GetMean();
979 rms = hisDelta->GetRMS();
987 "entries="<<entries<<
992 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1004 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1006 // Store object in the OCDB
1007 // By default the object is stored in the current directory
1008 // default comment consit of user name and the date
1010 TString ocdbStorage="";
1011 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1012 AliCDBMetaData *metaData= new AliCDBMetaData();
1013 metaData->SetObjectClassName("AliTPCCorrection");
1014 metaData->SetResponsible("Marian Ivanov");
1015 metaData->SetBeamPeriod(1);
1016 metaData->SetAliRootVersion("05-25-01"); //root version
1017 TString userName=gSystem->GetFromPipe("echo $USER");
1018 TString date=gSystem->GetFromPipe("date");
1020 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1021 if (comment) metaData->SetComment(comment);
1023 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1024 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1025 gStorage->Put(this, (*id1), metaData);