Added example file containing histograms of eta, pt and mult distributions.
[u/mrichter/AliRoot.git] / TPC / Base / 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
b4caed64 16// _________________________________________________________________
17//
18// Begin_Html
19// <h2> AliTPCCorrection class </h2>
20//
21// The AliTPCCorrection class provides a general framework to deal with space point distortions.
22// An correction class which inherits from here is for example AliTPCExBBShape or AliTPCExBTwist. <br>
23// General virtual functions are (for example) CorrectPoint(x,roc) where x is the vector of initial
24// positions in cartesian coordinates and roc represents the read-out chamber number according to
25// the offline numbering convention. The vector x is overwritten with the corrected coordinates. <br>
26// An alternative usage would be CorrectPoint(x,roc,dx), which leaves the vector x untouched, but
27// returns the distortions via the vector dx. <br>
28// This class is normally used via the general class AliTPCComposedCorrection.
29// <p>
30// Furthermore, the class contains basic geometrical descriptions like field cage radii
31// (fgkIFCRadius, fgkOFCRadius) and length (fgkTPCZ0) plus the voltages. Also, the definitions
32// of size and widths of the fulcrums building the grid of the final look-up table, which is
33// then interpolated, is defined in kNX and fgkXList).
34// <p>
35// All physics-model classes below are derived from this class in order to not duplicate code
36// and to allow a uniform treatment of all physics models.
37// <p>
38// <h3> Poisson solver </h3>
39// A numerical solver of the Poisson equation (relaxation technique) is implemented for 2-dimensional
40// geometries (r,z) as well as for 3-dimensional problems (r,$\phi$,z). The corresponding function
41// names are PoissonRelaxation?D. The relevant function arguments are the arrays of the boundary and
42// initial conditions (ArrayofArrayV, ArrayofChargeDensities) as well as the grid granularity which
43// is used during the calculation. These inputs can be chosen according to the needs of the physical
44// effect which is supposed to be simulated. In the 3D version, different symmetry conditions can be set
45// in order to reduce the calculation time (used in AliTPCFCVoltError3D).
46// <p>
47// <h3> Unified plotting functionality </h3>
48// Generic plot functions were implemented. They return a histogram pointer in the chosen plane of
49// the TPC drift volume with a selectable grid granularity and the magnitude of the correction vector.
50// For example, the function CreateHistoDZinXY(z,nx,ny) returns a 2-dimensional histogram which contains
51// the longitudinal corrections $dz$ in the (x,y)-plane at the given z position with the granularity of
52// nx and ny. The magnitude of the corrections is defined by the class from which this function is called.
53// In the same manner, standard plots for the (r,$\phi$)-plane and for the other corrections like $dr$ and $rd\phi$ are available
54// <p>
55// Note: This class is normally used via the class AliTPCComposedCorrection
56// End_Html
57//
6a1caa6b 58// Begin_Macro(source)
b4caed64 59// {
60// gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
6a1caa6b 61// TCanvas *c2 = new TCanvas("cAliTPCCorrection","cAliTPCCorrection",700,1050); c2->Divide(2,3);
b4caed64 62// AliTPCROCVoltError3D roc; // EXAMPLE PLOTS - SEE BELOW
63// roc.SetOmegaTauT1T2(0,1,1); // B=0
64// Float_t z0 = 1; // at +1 cm -> A side
65// c2->cd(1); roc.CreateHistoDRinXY(1.,300,300)->Draw("cont4z");
66// c2->cd(3);roc.CreateHistoDRPhiinXY(1.,300,300)->Draw("cont4z");
67// c2->cd(5);roc.CreateHistoDZinXY(1.,300,300)->Draw("cont4z");
68// Float_t phi0=0.5;
69// c2->cd(2);roc.CreateHistoDRinZR(phi0)->Draw("surf2");
70// c2->cd(4);roc.CreateHistoDRPhiinZR(phi0)->Draw("surf2");
71// c2->cd(6);roc.CreateHistoDZinZR(phi0)->Draw("surf2");
72// return c2;
73// }
74// End_Macro
75//
76// Begin_Html
77// <p>
78// Date: 27/04/2010 <br>
79// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
80// End_Html
81// _________________________________________________________________
82
83
be67055b 84#include "Riostream.h"
0116859c 85
86#include <TH2F.h>
87#include <TMath.h>
88#include <TROOT.h>
cf5b0aa0 89#include <TTreeStream.h>
ffab0c37 90#include <TTree.h>
91#include <TFile.h>
e527a1b9 92#include <TTimeStamp.h>
ffab0c37 93#include <AliCDBStorage.h>
94#include <AliCDBId.h>
95#include <AliCDBMetaData.h>
c9cbd2f2 96#include "TVectorD.h"
97#include "AliTPCParamSR.h"
7f4cb119 98
c9cbd2f2 99#include "AliTPCCorrection.h"
100#include "AliLog.h"
1b923461 101
1b923461 102#include "AliExternalTrackParam.h"
103#include "AliTrackPointArray.h"
104#include "TDatabasePDG.h"
105#include "AliTrackerBase.h"
106#include "AliTPCROC.h"
107#include "THnSparse.h"
108
c9cbd2f2 109#include "AliTPCLaserTrack.h"
110#include "AliESDVertex.h"
111#include "AliVertexerTracks.h"
112#include "TDatabasePDG.h"
113#include "TF1.h"
7f4cb119 114#include "TRandom.h"
c9cbd2f2 115
116#include "TDatabasePDG.h"
117
7f4cb119 118#include "AliTPCTransform.h"
119#include "AliTPCcalibDB.h"
120#include "AliTPCExB.h"
cf5b0aa0 121
c9cbd2f2 122#include "AliTPCRecoParam.h"
fdbbc146 123#include "TLinearFitter.h"
1b923461 124
0116859c 125
cf5b0aa0 126ClassImp(AliTPCCorrection)
127
f1817479 128
129TObjArray *AliTPCCorrection::fgVisualCorrection=0;
130// instance of correction for visualization
131
132
0116859c 133// FIXME: the following values should come from the database
c9cbd2f2 134const Double_t AliTPCCorrection::fgkTPCZ0 = 249.7; // nominal gating grid position
2b68ab9c 135const Double_t AliTPCCorrection::fgkIFCRadius= 83.5; // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
136// compare gkIFCRadius= 83.05: Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
c9cbd2f2 137const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
138const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
139const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
140const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
0116859c 141
c9cbd2f2 142const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
0116859c 143
c9cbd2f2 144const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
145const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
c9cbd2f2 146
0116859c 147
148AliTPCCorrection::AliTPCCorrection()
c9cbd2f2 149 : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 150{
151 //
152 // default constructor
153 //
f1817479 154 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
c9cbd2f2 155
35ae345f 156 InitLookUpfulcrums();
c9cbd2f2 157
0116859c 158}
159
160AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
c9cbd2f2 161: TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 162{
163 //
164 // default constructor, that set the name and title
165 //
f1817479 166 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
c9cbd2f2 167
35ae345f 168 InitLookUpfulcrums();
c9cbd2f2 169
0116859c 170}
171
172AliTPCCorrection::~AliTPCCorrection() {
173 //
174 // virtual destructor
175 //
176}
177
178void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
179 //
180 // Corrects the initial coordinates x (cartesian coordinates)
181 // according to the given effect (inherited classes)
182 // roc represents the TPC read out chamber (offline numbering convention)
183 //
184 Float_t dx[3];
185 GetCorrection(x,roc,dx);
186 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
187}
188
189void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
190 //
191 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
192 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
193 // roc represents the TPC read out chamber (offline numbering convention)
194 //
195 Float_t dx[3];
196 GetCorrection(x,roc,dx);
197 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
198}
199
200void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
201 //
202 // Distorts the initial coordinates x (cartesian coordinates)
203 // according to the given effect (inherited classes)
204 // roc represents the TPC read out chamber (offline numbering convention)
205 //
206 Float_t dx[3];
207 GetDistortion(x,roc,dx);
208 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
209}
210
46e89793 211void AliTPCCorrection::DistortPointLocal(Float_t x[],const Short_t roc) {
212 //
213 // Distorts the initial coordinates x (cartesian coordinates)
214 // according to the given effect (inherited classes)
215 // roc represents the TPC read out chamber (offline numbering convention)
216 //
217 Float_t gxyz[3]={0,0,0};
218 Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
219 Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
220 gxyz[0]= ca*x[0]+sa*x[1];
221 gxyz[1]= -sa*x[0]+ca*x[1];
222 gxyz[2]= x[2];
223 DistortPoint(gxyz,roc);
224 x[0]= ca*gxyz[0]-sa*gxyz[1];
225 x[1]= +sa*gxyz[0]+ca*gxyz[1];
226 x[2]= gxyz[2];
227}
228void AliTPCCorrection::CorrectPointLocal(Float_t x[],const Short_t roc) {
229 //
230 // Distorts the initial coordinates x (cartesian coordinates)
231 // according to the given effect (inherited classes)
232 // roc represents the TPC read out chamber (offline numbering convention)
233 //
234 Float_t gxyz[3]={0,0,0};
235 Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
236 Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
237 gxyz[0]= ca*x[0]+sa*x[1];
238 gxyz[1]= -sa*x[0]+ca*x[1];
239 gxyz[2]= x[2];
240 CorrectPoint(gxyz,roc);
241 x[0]= ca*gxyz[0]-sa*gxyz[1];
242 x[1]= sa*gxyz[0]+ca*gxyz[1];
243 x[2]= gxyz[2];
244}
245
0116859c 246void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
247 //
248 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
249 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
250 // roc represents the TPC read out chamber (offline numbering convention)
251 //
252 Float_t dx[3];
253 GetDistortion(x,roc,dx);
254 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
255}
256
257void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
258 //
259 // This function delivers the correction values dx in respect to the inital coordinates x
260 // roc represents the TPC read out chamber (offline numbering convention)
261 // Note: The dx is overwritten by the inherited effectice class ...
262 //
263 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
264}
265
266void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
267 //
268 // This function delivers the distortion values dx in respect to the inital coordinates x
269 // roc represents the TPC read out chamber (offline numbering convention)
270 //
271 GetCorrection(x,roc,dx);
272 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
273}
274
fdbbc146 275void AliTPCCorrection::GetCorrectionDz(const Float_t x[],const Short_t roc,Float_t dx[], Float_t delta) {
276 // author: marian.ivanov@cern.ch
277 //
278 // In this (virtual)function calculates the dx'/dz, dy'/dz and dz'/dz at given point (x,y,z)
279 // Generic implementation. Better precision can be acchieved knowing the internal structure
280 // of underlying trasnformation. Derived classes can reimplement it.
281 // To calculate correction is fitted in small neighberhood:
282 // (x+-delta,y+-delta,z+-delta) where delta is an argument
283 //
284 // Input parameters:
285 // x[] - space point corrdinate
286 // roc - readout chamber identifier (important e.g to do not miss the side of detector)
287 // delta - define the size of neighberhood
288 // Output parameter:
289 // dx[] - array {dx'/dz, dy'/dz , dz'/dz }
290
291 static TLinearFitter fitx(2,"pol1");
292 static TLinearFitter fity(2,"pol1");
293 static TLinearFitter fitz(2,"pol1");
294 fitx.ClearPoints();
295 fity.ClearPoints();
296 fitz.ClearPoints();
297 for (Int_t xdelta=-1; xdelta<=1; xdelta++)
298 for (Int_t ydelta=-1; ydelta<=1; ydelta++){
299 for (Int_t zdelta=-1; zdelta<=1; zdelta++){
300 Float_t xyz[3]={x[0]+xdelta*delta, x[1]+ydelta*delta, x[2]+zdelta*delta};
301 Float_t dxyz[3];
302 GetCorrection(xyz,roc,dxyz);
303 Double_t adelta=zdelta*delta;
304 fitx.AddPoint(&adelta, dxyz[0]);
305 fity.AddPoint(&adelta, dxyz[1]);
306 fitz.AddPoint(&adelta, dxyz[2]);
307 }
308 }
309 fitx.Eval();
310 fity.Eval();
311 fitz.Eval();
312 dx[0] = fitx.GetParameter(1);
313 dx[1] = fity.GetParameter(1);
314 dx[2] = fitz.GetParameter(1);
315}
316
317void AliTPCCorrection::GetCorrectionIntegralDz(const Float_t x[],const Short_t roc,Float_t dx[], Float_t delta){
318 //
319 // Integrate 3D distortion along drift lines
320 // To define the drift lines virtual function AliTPCCorrection::GetCorrectionDz is used
321 //
322 // Input parameters:
323 // x[] - space point corrdinate
324 // roc - readout chamber identifier (important e.g to do not miss the side of detector)
325 // delta - define the size of neighberhood
326 // Output parameter:
327 // dx[] - array { integral(dx'/dz), integral(dy'/dz) , integral(dz'/dz) }
328
329 Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
330 Double_t zdrift = TMath::Abs(x[2]-zroc);
331 Int_t nsteps = Int_t(zdrift/delta)+1;
332 //
333 //
334 Float_t xyz[3]={x[0],x[1],x[2]};
335 Float_t dxyz[3]={x[0],x[1],x[2]};
336 Float_t sign=((roc%36)<18) ? 1.:-1.;
337 Double_t sumdz=0;
338 for (Int_t i=0;i<nsteps; i++){
339 Float_t deltaZ=delta;
340 if (xyz[2]+deltaZ>fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
341 if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
342 deltaZ*=sign;
343 GetCorrectionDz(xyz,roc,dxyz,delta);
344 xyz[0]+=deltaZ*dxyz[0];
345 xyz[1]+=deltaZ*dxyz[1];
346 xyz[2]+=deltaZ; //
347 sumdz+=deltaZ*dxyz[2];
348 }
349 //
350 dx[0]=x[0]-xyz[0];
351 dx[1]=x[1]-xyz[1];
352 dx[2]= dxyz[2];
353
354}
355
05da1b4e 356void AliTPCCorrection::GetDistortionIntegralDz(const Float_t x[],const Short_t roc,Float_t dx[], Float_t delta){
357 //
358 // Integrate 3D distortion along drift lines
359 // To define the drift lines virtual function AliTPCCorrection::GetCorrectionDz is used
360 //
361 // Input parameters:
362 // x[] - space point corrdinate
363 // roc - readout chamber identifier (important e.g to do not miss the side of detector)
364 // delta - define the size of neighberhood
365 // Output parameter:
366 // dx[] - array { integral(dx'/dz), integral(dy'/dz) , integral(dz'/dz) }
367
368 Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
369 Double_t zdrift = TMath::Abs(x[2]-zroc);
370 Int_t nsteps = Int_t(zdrift/delta)+1;
371 //
372 //
373 Float_t xyz[3]={x[0],x[1],x[2]};
374 Float_t dxyz[3]={x[0],x[1],x[2]};
375 Float_t sign=((roc%36)<18) ? 1.:-1.;
376 Double_t sumdz=0;
377 for (Int_t i=0;i<nsteps; i++){
378 Float_t deltaZ=delta;
379 if (xyz[2]+deltaZ>fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
380 if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
381 deltaZ*=sign;
382 GetCorrectionDz(xyz,roc,dxyz,delta);
383 xyz[0]-=deltaZ*dxyz[0];
384 xyz[1]-=deltaZ*dxyz[1];
385 xyz[2]-=deltaZ; //
386 sumdz-=deltaZ*dxyz[2];
387 }
388 //
389 dx[0]=x[0]-xyz[0];
390 dx[1]=x[1]-xyz[1];
391 dx[2]= dxyz[2];
392
393}
fdbbc146 394
395
0116859c 396void AliTPCCorrection::Init() {
397 //
398 // Initialization funtion (not used at the moment)
399 //
400}
401
e527a1b9 402void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
403 //
404 // Update function
405 //
406}
407
0116859c 408void AliTPCCorrection::Print(Option_t* /*option*/) const {
409 //
410 // Print function to check which correction classes are used
411 // option=="d" prints details regarding the setted magnitude
412 // option=="a" prints the C0 and C1 coefficents for calibration purposes
413 //
414 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
415}
416
534fd34a 417void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
0116859c 418 //
419 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
420 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
421 // calibration run
422 //
534fd34a 423 fT1=t1;
424 fT2=t2;
425 //SetOmegaTauT1T2(omegaTau, t1, t2);
0116859c 426}
427
428TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
429 //
430 // Simple plot functionality.
431 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
432 // in respect to position z within the XY plane.
433 // The histogramm has nx times ny entries.
434 //
c9cbd2f2 435 AliTPCParam* tpcparam = new AliTPCParamSR;
436
0116859c 437 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
438 nx,-250.,250.,ny,-250.,250.);
439 Float_t x[3],dx[3];
440 x[2]=z;
441 Int_t roc=z>0.?0:18; // FIXME
442 for (Int_t iy=1;iy<=ny;++iy) {
443 x[1]=h->GetYaxis()->GetBinCenter(iy);
444 for (Int_t ix=1;ix<=nx;++ix) {
445 x[0]=h->GetXaxis()->GetBinCenter(ix);
446 GetCorrection(x,roc,dx);
447 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
c9cbd2f2 448 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
0116859c 449 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
450 h->SetBinContent(ix,iy,r1-r0);
451 }
452 else
453 h->SetBinContent(ix,iy,0.);
454 }
455 }
c9cbd2f2 456 delete tpcparam;
0116859c 457 return h;
458}
459
460TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
461 //
462 // Simple plot functionality.
463 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
464 // in respect to position z within the XY plane.
465 // The histogramm has nx times ny entries.
466 //
467
c9cbd2f2 468 AliTPCParam* tpcparam = new AliTPCParamSR;
469
0116859c 470 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
471 nx,-250.,250.,ny,-250.,250.);
472 Float_t x[3],dx[3];
473 x[2]=z;
474 Int_t roc=z>0.?0:18; // FIXME
475 for (Int_t iy=1;iy<=ny;++iy) {
476 x[1]=h->GetYaxis()->GetBinCenter(iy);
477 for (Int_t ix=1;ix<=nx;++ix) {
478 x[0]=h->GetXaxis()->GetBinCenter(ix);
479 GetCorrection(x,roc,dx);
480 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
c9cbd2f2 481 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
0116859c 482 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
483 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
484
485 Float_t dphi=phi1-phi0;
486 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
487 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
488
489 h->SetBinContent(ix,iy,r0*dphi);
490 }
491 else
492 h->SetBinContent(ix,iy,0.);
493 }
494 }
c9cbd2f2 495 delete tpcparam;
496 return h;
497}
498
499TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
500 //
501 // Simple plot functionality.
502 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
503 // in respect to position z within the XY plane.
504 // The histogramm has nx times ny entries.
505 //
506
507 AliTPCParam* tpcparam = new AliTPCParamSR;
508
509 TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
510 nx,-250.,250.,ny,-250.,250.);
511 Float_t x[3],dx[3];
512 x[2]=z;
513 Int_t roc=z>0.?0:18; // FIXME
514 for (Int_t iy=1;iy<=ny;++iy) {
515 x[1]=h->GetYaxis()->GetBinCenter(iy);
516 for (Int_t ix=1;ix<=nx;++ix) {
517 x[0]=h->GetXaxis()->GetBinCenter(ix);
518 GetCorrection(x,roc,dx);
519 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
520 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
521 h->SetBinContent(ix,iy,dx[2]);
522 }
523 else
524 h->SetBinContent(ix,iy,0.);
525 }
526 }
527 delete tpcparam;
0116859c 528 return h;
529}
530
531TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
532 //
533 // Simple plot functionality.
534 // Returns a 2d hisogram which represents the corrections in r direction (dr)
535 // in respect to angle phi within the ZR plane.
536 // The histogramm has nx times ny entries.
537 //
538 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
539 nz,-250.,250.,nr,85.,250.);
540 Float_t x[3],dx[3];
541 for (Int_t ir=1;ir<=nr;++ir) {
542 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
543 x[0]=radius*TMath::Cos(phi);
544 x[1]=radius*TMath::Sin(phi);
545 for (Int_t iz=1;iz<=nz;++iz) {
546 x[2]=h->GetXaxis()->GetBinCenter(iz);
547 Int_t roc=x[2]>0.?0:18; // FIXME
548 GetCorrection(x,roc,dx);
549 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
550 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
551 h->SetBinContent(iz,ir,r1-r0);
552 }
553 }
0116859c 554 return h;
555
556}
557
558TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
559 //
560 // Simple plot functionality.
561 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
562 // in respect to angle phi within the ZR plane.
563 // The histogramm has nx times ny entries.
564 //
565 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
566 nz,-250.,250.,nr,85.,250.);
567 Float_t x[3],dx[3];
568 for (Int_t iz=1;iz<=nz;++iz) {
569 x[2]=h->GetXaxis()->GetBinCenter(iz);
570 Int_t roc=x[2]>0.?0:18; // FIXME
571 for (Int_t ir=1;ir<=nr;++ir) {
572 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
573 x[0]=radius*TMath::Cos(phi);
574 x[1]=radius*TMath::Sin(phi);
575 GetCorrection(x,roc,dx);
576 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
577 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
578 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
579
580 Float_t dphi=phi1-phi0;
581 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
582 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
583
584 h->SetBinContent(iz,ir,r0*dphi);
585 }
586 }
587 return h;
588}
589
c9cbd2f2 590TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
591 //
592 // Simple plot functionality.
593 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
594 // in respect to angle phi within the ZR plane.
595 // The histogramm has nx times ny entries.
596 //
597 TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
598 nz,-250.,250.,nr,85.,250.);
599 Float_t x[3],dx[3];
600 for (Int_t ir=1;ir<=nr;++ir) {
601 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
602 x[0]=radius*TMath::Cos(phi);
603 x[1]=radius*TMath::Sin(phi);
604 for (Int_t iz=1;iz<=nz;++iz) {
605 x[2]=h->GetXaxis()->GetBinCenter(iz);
606 Int_t roc=x[2]>0.?0:18; // FIXME
607 GetCorrection(x,roc,dx);
608 h->SetBinContent(iz,ir,dx[2]);
609 }
610 }
611 return h;
612
613}
614
615
0116859c 616TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
617 const char *xlabel,const char *ylabel,const char *zlabel,
618 Int_t nbinsx,Double_t xlow,Double_t xup,
619 Int_t nbinsy,Double_t ylow,Double_t yup) {
620 //
621 // Helper function to create a 2d histogramm of given size
622 //
623
624 TString hname=name;
625 Int_t i=0;
626 if (gDirectory) {
627 while (gDirectory->FindObject(hname.Data())) {
628 hname =name;
629 hname+="_";
630 hname+=i;
631 ++i;
632 }
633 }
634 TH2F *h=new TH2F(hname.Data(),title,
635 nbinsx,xlow,xup,
636 nbinsy,ylow,yup);
637 h->GetXaxis()->SetTitle(xlabel);
638 h->GetYaxis()->SetTitle(ylabel);
639 h->GetZaxis()->SetTitle(zlabel);
640 h->SetStats(0);
641 return h;
642}
643
0116859c 644// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
645
646void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
b1f0a2a5 647 const Double_t er[kNZ][kNR], Double_t &erValue ) {
0116859c 648 //
649 // Interpolate table - 2D interpolation
650 //
25732bff 651 Double_t saveEr[5] = {0,0,0,0,0};
0116859c 652
653 Search( kNZ, fgkZList, z, fJLow ) ;
654 Search( kNR, fgkRList, r, fKLow ) ;
655 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
656 if ( fKLow < 0 ) fKLow = 0 ;
657 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
658 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
659
660 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
b1f0a2a5 661 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
0116859c 662 }
b1f0a2a5 663 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
0116859c 664
665}
666
c9cbd2f2 667void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
668 const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
669 Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
670 //
671 // Interpolate table - 3D interpolation
672 //
673
25732bff 674 Double_t saveEr[5]= {0,0,0,0,0};
675 Double_t savedEr[5]= {0,0,0,0,0} ;
676
677 Double_t saveEphi[5]= {0,0,0,0,0};
678 Double_t savedEphi[5]= {0,0,0,0,0} ;
679
680 Double_t saveEz[5]= {0,0,0,0,0};
681 Double_t savedEz[5]= {0,0,0,0,0} ;
c9cbd2f2 682
683 Search( kNZ, fgkZList, z, fILow ) ;
684 Search( kNPhi, fgkPhiList, z, fJLow ) ;
685 Search( kNR, fgkRList, r, fKLow ) ;
686
687 if ( fILow < 0 ) fILow = 0 ; // check if out of range
688 if ( fJLow < 0 ) fJLow = 0 ;
689 if ( fKLow < 0 ) fKLow = 0 ;
690
691 if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
692 if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
693 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
694
695 for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
696 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
697 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
698 saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
699 saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
700 }
701 savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
702 savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
703 savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
704 }
705 erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
706 ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
707 ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
708
709}
710
711Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
712 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
713 const TMatrixD &array ) {
714 //
715 // Interpolate table (TMatrix format) - 2D interpolation
716 //
717
718 static Int_t jlow = 0, klow = 0 ;
25732bff 719 Double_t saveArray[5] = {0,0,0,0,0} ;
c9cbd2f2 720
721 Search( nx, xv, x, jlow ) ;
722 Search( ny, yv, y, klow ) ;
723 if ( jlow < 0 ) jlow = 0 ; // check if out of range
724 if ( klow < 0 ) klow = 0 ;
725 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
726 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
727
728 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
729 {
730 Double_t *ajkl = &((TMatrixD&)array)(j,klow);
731 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
732 }
733
734 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
735
736}
737
738Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
739 const Int_t nx, const Int_t ny, const Int_t nz,
740 const Double_t xv[], const Double_t yv[], const Double_t zv[],
741 TMatrixD **arrayofArrays ) {
742 //
743 // Interpolate table (TMatrix format) - 3D interpolation
744 //
745
746 static Int_t ilow = 0, jlow = 0, klow = 0 ;
25732bff 747 Double_t saveArray[5]= {0,0,0,0,0};
748 Double_t savedArray[5]= {0,0,0,0,0} ;
c9cbd2f2 749
750 Search( nx, xv, x, ilow ) ;
751 Search( ny, yv, y, jlow ) ;
752 Search( nz, zv, z, klow ) ;
753
754 if ( ilow < 0 ) ilow = 0 ; // check if out of range
755 if ( jlow < 0 ) jlow = 0 ;
756 if ( klow < 0 ) klow = 0 ;
757
758 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
759 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
760 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
761
762 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
763 {
764 TMatrixD &table = *arrayofArrays[k] ;
765 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
766 {
767 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
768 }
769 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
770 }
771 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
772
773}
774
0116859c 775Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 776 const Int_t order, const Double_t x ) {
0116859c 777 //
778 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
779 //
780
781 Double_t y ;
782 if ( order == 2 ) { // Quadratic Interpolation = 2
783 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
784 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
785 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
786 } else { // Linear Interpolation = 1
787 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
788 }
789
790 return (y);
791
792}
793
2bf29b72 794Float_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
795 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
796 const TMatrixF &array ) {
797 //
798 // Interpolate table (TMatrix format) - 2D interpolation
799 // Float version (in order to decrease the OCDB size)
800 //
801
802 static Int_t jlow = 0, klow = 0 ;
803 Float_t saveArray[5] = {0.,0.,0.,0.,0.} ;
804
805 Search( nx, xv, x, jlow ) ;
806 Search( ny, yv, y, klow ) ;
807 if ( jlow < 0 ) jlow = 0 ; // check if out of range
808 if ( klow < 0 ) klow = 0 ;
809 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
810 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
811
812 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
813 {
814 Float_t *ajkl = &((TMatrixF&)array)(j,klow);
815 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
816 }
817
818 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
819
820}
821
822Float_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
823 const Int_t nx, const Int_t ny, const Int_t nz,
824 const Double_t xv[], const Double_t yv[], const Double_t zv[],
825 TMatrixF **arrayofArrays ) {
826 //
827 // Interpolate table (TMatrix format) - 3D interpolation
828 // Float version (in order to decrease the OCDB size)
829 //
830
831 static Int_t ilow = 0, jlow = 0, klow = 0 ;
832 Float_t saveArray[5]= {0.,0.,0.,0.,0.};
833 Float_t savedArray[5]= {0.,0.,0.,0.,0.} ;
834
835 Search( nx, xv, x, ilow ) ;
836 Search( ny, yv, y, jlow ) ;
837 Search( nz, zv, z, klow ) ;
838
839 if ( ilow < 0 ) ilow = 0 ; // check if out of range
840 if ( jlow < 0 ) jlow = 0 ;
841 if ( klow < 0 ) klow = 0 ;
842
843 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
844 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
845 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
846
847 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
848 {
849 TMatrixF &table = *arrayofArrays[k] ;
850 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
851 {
852 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
853 }
854 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
855 }
856 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
857
858}
859Float_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Float_t yArray[],
860 const Int_t order, const Double_t x ) {
861 //
862 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
863 // Float version (in order to decrease the OCDB size)
864 //
865
866 Float_t y ;
867 if ( order == 2 ) { // Quadratic Interpolation = 2
868 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
869 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
870 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
871 } else { // Linear Interpolation = 1
872 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
873 }
874
875 return (y);
876
877}
878
879
0116859c 880
b1f0a2a5 881void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 882 //
883 // Search an ordered table by starting at the most recently used point
884 //
885
886 Long_t middle, high ;
887 Int_t ascend = 0, increment = 1 ;
888
889 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
890
891 if ( low < 0 || low > n-1 ) {
892 low = -1 ; high = n ;
893 } else { // Ordered Search phase
894 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
895 if ( low == n-1 ) return ;
896 high = low + 1 ;
897 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
898 low = high ;
899 increment *= 2 ;
900 high = low + increment ;
901 if ( high > n-1 ) { high = n ; break ; }
902 }
903 } else {
904 if ( low == 0 ) { low = -1 ; return ; }
905 high = low - 1 ;
906 while ( (Int_t)( x < xArray[low] ) == ascend ) {
907 high = low ;
908 increment *= 2 ;
909 if ( increment >= high ) { low = -1 ; break ; }
910 else low = high - increment ;
911 }
912 }
913 }
914
915 while ( (high-low) != 1 ) { // Binary Search Phase
916 middle = ( high + low ) / 2 ;
917 if ( (Int_t)( x >= xArray[middle] ) == ascend )
918 low = middle ;
919 else
920 high = middle ;
921 }
922
923 if ( x == xArray[n-1] ) low = n-2 ;
924 if ( x == xArray[0] ) low = 0 ;
925
926}
927
35ae345f 928void AliTPCCorrection::InitLookUpfulcrums() {
929 //
930 // Initialization of interpolation points - for main look up table
931 // (course grid in the middle, fine grid on the borders)
932 //
933
934 AliTPCROC * roc = AliTPCROC::Instance();
935 const Double_t rLow = TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
936
937 // fulcrums in R
938 fgkRList[0] = rLow;
939 for (Int_t i = 1; i<kNR; i++) {
940 fgkRList[i] = fgkRList[i-1] + 3.5; // 3.5 cm spacing
941 if (fgkRList[i]<90 ||fgkRList[i]>245)
942 fgkRList[i] = fgkRList[i-1] + 0.5; // 0.5 cm spacing
943 else if (fgkRList[i]<100 || fgkRList[i]>235)
944 fgkRList[i] = fgkRList[i-1] + 1.5; // 1.5 cm spacing
945 else if (fgkRList[i]<120 || fgkRList[i]>225)
946 fgkRList[i] = fgkRList[i-1] + 2.5; // 2.5 cm spacing
947 }
948
949 // fulcrums in Z
950 fgkZList[0] = -249.5;
951 fgkZList[kNZ-1] = 249.5;
952 for (Int_t j = 1; j<kNZ/2; j++) {
953 fgkZList[j] = fgkZList[j-1];
954 if (TMath::Abs(fgkZList[j])< 0.15)
955 fgkZList[j] = fgkZList[j-1] + 0.09; // 0.09 cm spacing
956 else if(TMath::Abs(fgkZList[j])< 0.6)
957 fgkZList[j] = fgkZList[j-1] + 0.4; // 0.4 cm spacing
958 else if (TMath::Abs(fgkZList[j])< 2.5 || TMath::Abs(fgkZList[j])>248)
959 fgkZList[j] = fgkZList[j-1] + 0.5; // 0.5 cm spacing
960 else if (TMath::Abs(fgkZList[j])<10 || TMath::Abs(fgkZList[j])>235)
961 fgkZList[j] = fgkZList[j-1] + 1.5; // 1.5 cm spacing
962 else if (TMath::Abs(fgkZList[j])<25 || TMath::Abs(fgkZList[j])>225)
963 fgkZList[j] = fgkZList[j-1] + 2.5; // 2.5 cm spacing
964 else
965 fgkZList[j] = fgkZList[j-1] + 4; // 4 cm spacing
966
967 fgkZList[kNZ-j-1] = -fgkZList[j];
968 }
969
970 // fulcrums in phi
971 for (Int_t k = 0; k<kNPhi; k++)
972 fgkPhiList[k] = TMath::TwoPi()*k/(kNPhi-1);
973
974
975}
976
977
c9cbd2f2 978void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
979 TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
980 const Int_t rows, const Int_t columns, const Int_t iterations,
981 const Bool_t rocDisplacement ) {
1b923461 982 //
983 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
984 //
985 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
986 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
987 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
988 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
989 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
990 // random numbers back as a solution.
991 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
992 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
993 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
994 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
995 // number of points and solves the problem again. This happens several times until the maximum number
996 // of points has been included in the array.
997 //
998 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
999 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
1000 //
c9cbd2f2 1001 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
1002 //
1b923461 1003 // Original code by Jim Thomas (STAR TPC Collaboration)
1004 //
1005
1006 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
1007
1008 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
1009 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
1010 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
1011
1012 TMatrixD arrayEr(rows,columns) ;
1013 TMatrixD arrayEz(rows,columns) ;
1014
1015 //Check that number of rows and columns is suitable for a binary expansion
1016
1017 if ( !IsPowerOfTwo(rows-1) ) {
1018 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
1019 return;
1020 }
1021 if ( !IsPowerOfTwo(columns-1) ) {
1022 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
1023 return;
1024 }
1025
1026 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
1027 // Allow for different size grid spacing in R and Z directions
1028 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
1029
1030 Int_t iOne = (rows-1)/4 ;
1031 Int_t jOne = (columns-1)/4 ;
1032 // Solve for N in 2**N, add one.
1033 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
1034
1035 for ( Int_t count = 0 ; count < loops ; count++ ) {
1036 // Loop while the matrix expands & the resolution increases.
1037
1038 Float_t tempGridSizeR = gridSizeR * iOne ;
1039 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
1040 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
1041
1042 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1043 std::vector<float> coef1(rows) ;
1044 std::vector<float> coef2(rows) ;
1045
1046 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
1047 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1048 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
1049 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
1050 }
1051
1052 TMatrixD sumChargeDensity(rows,columns) ;
1053
1054 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1055 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
1056 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1057 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1058 else {
1059 // Add up all enclosed charge density contributions within 1/2 unit in all directions
1060 Float_t weight = 0.0 ;
1061 Float_t sum = 0.0 ;
1062 sumChargeDensity(i,j) = 0.0 ;
1063 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
1064 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
1065 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
1066 else
1067 weight = 1.0 ;
1068 // Note that this is cylindrical geometry
1069 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1070 sum += weight*radius ;
1071 }
1072 }
1073 sumChargeDensity(i,j) /= sum ;
1074 }
1075 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
1076 }
1077 }
1078
1079 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1080 // Solve Poisson's Equation
1081 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
1082 // as interations increase.
1083 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1084 Float_t overRelaxM1 = overRelax - 1.0 ;
1085 Float_t overRelaxtempFourth, overRelaxcoef5 ;
1086 overRelaxtempFourth = overRelax * tempFourth ;
1087 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
1088
1089 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1090 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1091
1092 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
1093 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
1094 - overRelaxcoef5 * arrayV(i,j)
1095 + coef1[i] * arrayV(i+iOne,j)
1096 + sumChargeDensity(i,j)
1097 ) * overRelaxtempFourth;
1098 }
1099 }
1100
1101 if ( k == iterations ) {
1102 // After full solution is achieved, copy low resolution solution into higher res array
1103 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1104 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1105
1106 if ( iOne > 1 ) {
1107 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
1108 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
1109 }
1110 if ( jOne > 1 ) {
1111 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
1112 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
1113 }
1114 if ( iOne > 1 && jOne > 1 ) {
1115 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
1116 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
1117 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
1118 // Note that this leaves a point at the upper left and lower right corners uninitialized.
1119 // -> Not a big deal.
1120 }
1121
1122 }
1123 }
1124 }
1125
1126 }
1127
1128 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
1129 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
1130
c9cbd2f2 1131 sumChargeDensity.Clear();
1b923461 1132 }
1133
1134 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
1135 for ( Int_t j = 0 ; j < columns ; j++ ) {
1136 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1137 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1138 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1139 }
1140
1141 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1142 for ( Int_t i = 0 ; i < rows ; i++) {
1143 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1144 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1145 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1146 }
1147
1148 for ( Int_t i = 0 ; i < rows ; i++) {
1149 // Note: go back and compare to old version of this code. See notes below.
1150 // JT Test ... attempt to divide by real Ez not Ez to first order
1151 for ( Int_t j = 0 ; j < columns ; j++ ) {
1152 arrayEz(i,j) += ezField;
1153 // This adds back the overall Z gradient of the field (main E field component)
1154 }
1155 // Warning: (-=) assumes you are using an error potetial without the overall Field included
1156 }
1157
1158 // Integrate Er/Ez from Z to zero
1159 for ( Int_t j = 0 ; j < columns ; j++ ) {
1160 for ( Int_t i = 0 ; i < rows ; i++ ) {
c9cbd2f2 1161
1b923461 1162 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1163 arrayErOverEz(i,j) = 0.0 ;
c9cbd2f2 1164 arrayDeltaEz(i,j) = 0.0 ;
1165
1b923461 1166 for ( Int_t k = j ; k < columns ; k++ ) {
1167 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
c9cbd2f2 1168 arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
1b923461 1169 if ( index != 4 ) index = 4; else index = 2 ;
1170 }
c9cbd2f2 1171 if ( index == 4 ) {
1172 arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
1173 arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
1174 }
1175 if ( index == 2 ) {
1176 arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1177 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
1178 arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
1179 -2.5*(arrayEz(i,columns-1)-ezField));
1180 }
1181 if ( j == columns-2 ) {
1182 arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1183 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
1184 arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
1185 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
1186 }
1187 if ( j == columns-1 ) {
1188 arrayErOverEz(i,j) = 0.0 ;
1189 arrayDeltaEz(i,j) = 0.0 ;
1190 }
1b923461 1191 }
1192 }
1193
c9cbd2f2 1194 // calculate z distortion from the integrated Delta Ez residuals
1195 // and include the aquivalence (Volt to cm) of the ROC shift !!
1196
1197 for ( Int_t j = 0 ; j < columns ; j++ ) {
1198 for ( Int_t i = 0 ; i < rows ; i++ ) {
1199
1200 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1201 arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
1202
1203 // ROC Potential in cm aquivalent
1204 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1205 if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
1206
1207 }
1208 }
1209
1210 arrayEr.Clear();
1211 arrayEz.Clear();
1212
1b923461 1213}
1214
c9cbd2f2 1215void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
1216 TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
1217 const Int_t rows, const Int_t columns, const Int_t phislices,
1218 const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
1219 Bool_t rocDisplacement ) {
1220 //
1221 // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
1222 //
1223 // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
1224 // The number of rows and COLUMNS can be different.
1225 //
1226 // ROWS == 2**M + 1
1227 // COLUMNS == 2**N + 1
1228 // PHISLICES == Arbitrary but greater than 3
1229 //
1230 // DeltaPhi in Radians
1231 //
1232 // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
1233 // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
1234 //
1235 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
1236
1237 const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
1238
1239 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
1240 const Float_t gridSizePhi = deltaphi ;
1241 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
1242 const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
1243 const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
1244
1245 TMatrixD arrayE(rows,columns) ;
1246
1247 // Check that the number of rows and columns is suitable for a binary expansion
1248 if ( !IsPowerOfTwo((rows-1)) ) {
1249 AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
1250 return; }
1251 if ( !IsPowerOfTwo((columns-1)) ) {
1252 AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
1253 return; }
1254 if ( phislices <= 3 ) {
1255 AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
1256 return; }
1257 if ( phislices > 1000 ) {
1258 AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
1259 return; }
1260
1261 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
1262 // Allow for different size grid spacing in R and Z directions
1263 // Use a binary expansion of the matrix to speed up the solution of the problem
1264
1265 Int_t loops, mplus, mminus, signplus, signminus ;
1266 Int_t ione = (rows-1)/4 ;
1267 Int_t jone = (columns-1)/4 ;
1268 loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
1269 loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
1270
1271 TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
1272
1273 for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
1274
1275 for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
1276
1277 Float_t tempgridSizeR = gridSizeR * ione ;
1278 Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
1279 Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
1280
1281 std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1282 std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1283 std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1284 std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1285
1286 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1287 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1288 coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1289 coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1290 coef3[i] = tempratioPhi/(radius*radius);
1291 coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1292 }
1293
1294 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1295 TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
1296 TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1297 for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1298 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1299 for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1300 if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1301 else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
1302 Float_t weight = 0.0 ;
1303 Float_t sum = 0.0 ;
1304 sumChargeDensity(i,j) = 0.0 ;
1305 for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1306 for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1307 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1308 else
1309 weight = 1.0 ;
1310 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1311 sum += weight*radius ;
1312 }
1313 }
1314 sumChargeDensity(i,j) /= sum ;
1315 }
1316 sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1317 }
1318 }
1319 }
1320
1321 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1322
1323 // over-relaxation index, >= 1 but < 2
1324 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1325 Float_t overRelaxM1 = overRelax - 1.0 ;
1326
1327 std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
1328 std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
1329
1330 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1331 overRelaxcoef4[i] = overRelax * coef4[i] ;
1332 overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1333 }
1334
1335 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1336
1337 mplus = m + 1; signplus = 1 ;
1338 mminus = m - 1 ; signminus = 1 ;
1339 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1340 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1341 if ( mminus < 0 ) mminus = 1 ;
1342 }
1343 else if (symmetry==-1) { // Anti-symmetry in phi
1344 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1345 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1346 }
1347 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1348 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1349 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1350 }
1351 TMatrixD& arrayV = *arrayofArrayV[m] ;
1352 TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
1353 TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
1354 TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
1355
1356 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1357 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1358
1359 arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
1360 + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1361 - overRelaxcoef5[i] * arrayV(i,j)
1362 + coef1[i] * arrayV(i+ione,j)
1363 + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1364 + sumChargeDensity(i,j)
1365 ) * overRelaxcoef4[i] ;
1366 // Note: over-relax the solution at each step. This speeds up the convergance.
1367
1368 }
1369 }
1370
1371 if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
1372 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1373 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1374
1375 if ( ione > 1 ) {
1376 arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
1377 if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
1378 }
1379 if ( jone > 1 ) {
1380 arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
1381 if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
1382 }
1383 if ( ione > 1 && jone > 1 ) {
1384 arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1385 if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1386 if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1387 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1388 }
1389 }
1390 }
1391 }
1392
1393 }
1394 }
1395
1396 ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1397 jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1398
1399 }
1400
1401 //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1402 //Integrate E(r)/E(z) from point of origin to pad plane
1403
1404 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1405 TMatrixD& arrayV = *arrayofArrayV[m] ;
1406 TMatrixD& eroverEz = *arrayofEroverEz[m] ;
1407
1408 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1409
1410 // Differentiate in R
1411 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1412 arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1413 arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1414 // Integrate over Z
1415 for ( Int_t i = 0 ; i < rows ; i++ ) {
1416 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1417 eroverEz(i,j) = 0.0 ;
1418 for ( Int_t k = j ; k < columns ; k++ ) {
1419
1420 eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1421 if ( index != 4 ) index = 4; else index = 2 ;
1422 }
1423 if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1424 if ( index == 2 ) eroverEz(i,j) +=
1425 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1426 if ( j == columns-2 ) eroverEz(i,j) =
1427 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1428 if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
1429 }
1430 }
1431 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1432 // eroverEz.Draw("surf") ; } // JT test
1433 }
1434
1435 //Differentiate V(r) and solve for E(phi)
1436 //Integrate E(phi)/E(z) from point of origin to pad plane
1437
1438 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1439
1440 mplus = m + 1; signplus = 1 ;
1441 mminus = m - 1 ; signminus = 1 ;
1442 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1443 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1444 if ( mminus < 0 ) mminus = 1 ;
1445 }
1446 else if (symmetry==-1) { // Anti-symmetry in phi
1447 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1448 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1449 }
1450 else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1451 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1452 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1453 }
1454 TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
1455 TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
1456 TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
1457 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1458 // Differentiate in Phi
1459 for ( Int_t i = 0 ; i < rows ; i++ ) {
1460 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1461 arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1462 }
1463 // Integrate over Z
1464 for ( Int_t i = 0 ; i < rows ; i++ ) {
1465 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1466 ePhioverEz(i,j) = 0.0 ;
1467 for ( Int_t k = j ; k < columns ; k++ ) {
1468
1469 ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1470 if ( index != 4 ) index = 4; else index = 2 ;
1471 }
1472 if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1473 if ( index == 2 ) ePhioverEz(i,j) +=
1474 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1475 if ( j == columns-2 ) ePhioverEz(i,j) =
1476 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1477 if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
1478 }
1479 }
1480 // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
1481 // arrayE.Draw("surf") ; } // JT test
1482 }
1483
1484
1485 // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1486 // Integrate (E(z)-Ezstd) from point of origin to pad plane
1487
1488 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1489 TMatrixD& arrayV = *arrayofArrayV[m] ;
1490 TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
1491
1492 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1493 for ( Int_t i = 0 ; i < rows ; i++) {
1494 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1495 arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1496 arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1497 }
1498
1499 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1500 // Integrate over Z
1501 for ( Int_t i = 0 ; i < rows ; i++ ) {
1502 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1503 deltaEz(i,j) = 0.0 ;
1504 for ( Int_t k = j ; k < columns ; k++ ) {
1505 deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
1506 if ( index != 4 ) index = 4; else index = 2 ;
1507 }
1508 if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1509 if ( index == 2 ) deltaEz(i,j) +=
1510 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1511 if ( j == columns-2 ) deltaEz(i,j) =
1512 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1513 if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
1514 }
1515 }
1516 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1517 // eroverEz.Draw("surf") ; } // JT test
1518
1519 // calculate z distortion from the integrated Delta Ez residuals
1520 // and include the aquivalence (Volt to cm) of the ROC shift !!
1521
1522 for ( Int_t j = 0 ; j < columns ; j++ ) {
1523 for ( Int_t i = 0 ; i < rows ; i++ ) {
1524
1525 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1526 deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1527
1528 // ROC Potential in cm aquivalent
1529 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1530 if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
1531
1532 }
1533 }
1534
1535 } // end loop over phi
1536
1537
1538
1539 for ( Int_t k = 0 ; k < phislices ; k++ )
1540 {
1541 arrayofSumChargeDensities[k]->Delete() ;
1542 }
1543
1544
1545
1546 arrayE.Clear();
1547}
1b923461 1548
1549
710bda39 1550Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1b923461 1551 //
1552 // Helperfunction: Check if integer is a power of 2
1553 //
1554 Int_t j = 0;
1555 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1556 if ( j == 1 ) return(1) ; // True
1557 return(0) ; // False
1558}
1559
cf5b0aa0 1560
b1f0a2a5 1561AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 1562 //
1563 // Fit the track parameters - without and with distortion
1564 // 1. Space points in the TPC are simulated along the trajectory
1565 // 2. Space points distorted
1566 // 3. Fits the non distorted and distroted track to the reference plane at refX
1567 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
1568 //
1569 // trackIn - input track parameters
1570 // refX - reference X to fit the track
1571 // dir - direction - out=1 or in=-1
1572 // pcstream - debug streamer to check the results
1573 //
cad404e1 1574 // see AliExternalTrackParam.h documentation:
1575 // track1.fP[0] - local y (rphi)
1576 // track1.fP[1] - z
1577 // track1.fP[2] - sinus of local inclination angle
1578 // track1.fP[3] - tangent of deep angle
1579 // track1.fP[4] - 1/pt
1b923461 1580
cf5b0aa0 1581 AliTPCROC * roc = AliTPCROC::Instance();
1582 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1583 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
1584 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
cf5b0aa0 1585 const Double_t kMaxSnp = 0.85;
1586 const Double_t kSigmaY=0.1;
1587 const Double_t kSigmaZ=0.1;
ca58ed4e 1588 const Double_t kMaxR=500;
1589 const Double_t kMaxZ=500;
46e89793 1590
cfe2c39a 1591 const Double_t kMaxZ0=220;
1592 const Double_t kZcut=3;
cf5b0aa0 1593 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
ca58ed4e 1594 Int_t npoints1=0;
1595 Int_t npoints2=0;
cf5b0aa0 1596
be67055b 1597 AliExternalTrackParam track(trackIn); //
cf5b0aa0 1598 // generate points
1599 AliTrackPointArray pointArray0(npoints0);
1600 AliTrackPointArray pointArray1(npoints0);
1601 Double_t xyz[3];
cfe2c39a 1602 if (!AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,5,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1603 //
1604 // simulate the track
1605 Int_t npoints=0;
1606 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
1607 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
cfe2c39a 1608 if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1609 track.GetXYZ(xyz);
cfe2c39a 1610 xyz[0]+=gRandom->Gaus(0,0.000005);
1611 xyz[1]+=gRandom->Gaus(0,0.000005);
1612 xyz[2]+=gRandom->Gaus(0,0.000005);
1613 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
46e89793 1614 if (TMath::Abs(track.GetX())<kRTPC0) continue;
1615 if (TMath::Abs(track.GetX())>kRTPC1) continue;
cf5b0aa0 1616 AliTrackPoint pIn0; // space point
1617 AliTrackPoint pIn1;
ffab0c37 1618 Int_t sector= (xyz[2]>0)? 0:18;
cf5b0aa0 1619 pointArray0.GetPoint(pIn0,npoints);
1620 pointArray1.GetPoint(pIn1,npoints);
1621 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1622 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
ffab0c37 1623 DistortPoint(distPoint, sector);
cf5b0aa0 1624 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1625 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1626 //
1627 track.Rotate(alpha);
1628 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
1629 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
1630 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1631 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1632 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
1633 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
1634 pointArray0.AddPoint(npoints, &pIn0);
1635 pointArray1.AddPoint(npoints, &pIn1);
1636 npoints++;
1637 if (npoints>=npoints0) break;
1638 }
cfe2c39a 1639 if (npoints<npoints0/4.) return 0;
cf5b0aa0 1640 //
1641 // refit track
1642 //
1643 AliExternalTrackParam *track0=0;
1644 AliExternalTrackParam *track1=0;
1645 AliTrackPoint point1,point2,point3;
1646 if (dir==1) { //make seed inner
1647 pointArray0.GetPoint(point1,1);
cfe2c39a 1648 pointArray0.GetPoint(point2,11);
1649 pointArray0.GetPoint(point3,21);
cf5b0aa0 1650 }
1651 if (dir==-1){ //make seed outer
cfe2c39a 1652 pointArray0.GetPoint(point1,npoints-21);
1653 pointArray0.GetPoint(point2,npoints-11);
cf5b0aa0 1654 pointArray0.GetPoint(point3,npoints-1);
46e89793 1655 }
1656 if ((TMath::Abs(point1.GetX()-point3.GetX())+TMath::Abs(point1.GetY()-point3.GetY()))<10){
1657 printf("fit points not properly initialized\n");
1658 return 0;
1659 }
cf5b0aa0 1660 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1661 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
cfe2c39a 1662 track0->ResetCovariance(10);
1663 track1->ResetCovariance(10);
1664 if (TMath::Abs(AliTrackerBase::GetBz())<0.01){
1665 ((Double_t*)track0->GetParameter())[4]= trackIn.GetParameter()[4];
1666 ((Double_t*)track1->GetParameter())[4]= trackIn.GetParameter()[4];
1667 }
cf5b0aa0 1668 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
8b63d99c 1669 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
cf5b0aa0 1670 //
1671 AliTrackPoint pIn0;
1672 AliTrackPoint pIn1;
1673 pointArray0.GetPoint(pIn0,ipoint);
1674 pointArray1.GetPoint(pIn1,ipoint);
1675 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
1676 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
46e89793 1677 if (TMath::Abs(prot0.GetX())<kRTPC0) continue;
1678 if (TMath::Abs(prot0.GetX())>kRTPC1) continue;
cf5b0aa0 1679 //
cfe2c39a 1680 if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
1681 if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
ca58ed4e 1682 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1683 if (TMath::Abs(track0->GetX())>kMaxR) break;
1684 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1685 if (TMath::Abs(track1->GetX())>kMaxR) break;
cfe2c39a 1686 if (dir>0 && track1->GetX()>refX) continue;
1687 if (dir<0 && track1->GetX()<refX) continue;
1688 if (TMath::Abs(track1->GetZ())<kZcut)continue;
8b63d99c 1689 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
cf5b0aa0 1690 //
1691 Double_t pointPos[2]={0,0};
1692 Double_t pointCov[3]={0,0,0};
1693 pointPos[0]=prot0.GetY();//local y
1694 pointPos[1]=prot0.GetZ();//local z
1695 pointCov[0]=prot0.GetCov()[3];//simay^2
1696 pointCov[1]=prot0.GetCov()[4];//sigmayz
1697 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
ca58ed4e 1698 if (!track0->Update(pointPos,pointCov)) break;
cf5b0aa0 1699 //
8b63d99c 1700 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
1701 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
1702 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
1703
0b736a46 1704 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1705 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
cf5b0aa0 1706 pointCov[0]=prot1.GetCov()[3];//simay^2
1707 pointCov[1]=prot1.GetCov()[4];//sigmayz
1708 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
ca58ed4e 1709 if (!track1->Update(pointPos,pointCov)) break;
1710 npoints1++;
1711 npoints2++;
cf5b0aa0 1712 }
cfe2c39a 1713 if (npoints2<npoints/4.) return 0;
1714 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,5.,kTRUE,kMaxSnp);
1715 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,1.,kTRUE,kMaxSnp);
cf5b0aa0 1716 track1->Rotate(track0->GetAlpha());
cfe2c39a 1717 AliTrackerBase::PropagateTrackTo(track1,track0->GetX(),kMass,5.,kFALSE,kMaxSnp);
cf5b0aa0 1718
cad404e1 1719 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 1720 "point0.="<<&pointArray0<< // points
1721 "point1.="<<&pointArray1<< // distorted points
1722 "trackIn.="<<&track<< // original track
1723 "track0.="<<track0<< // fitted track
1724 "track1.="<<track1<< // fitted distorted track
1725 "\n";
be67055b 1726 new(&trackIn) AliExternalTrackParam(*track0);
cf5b0aa0 1727 delete track0;
1728 return track1;
1729}
1730
1731
ffab0c37 1732
1733
1734
1735TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1736 //
1737 // create the distortion tree on a mesh with granularity given by step
1738 // return the tree with distortions at given position
1739 // Map is created on the mesh with given step size
1740 //
1741 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1742 Float_t xyz[3];
1743 for (Double_t x= -250; x<250; x+=step){
1744 for (Double_t y= -250; y<250; y+=step){
1745 Double_t r = TMath::Sqrt(x*x+y*y);
1746 if (r<80) continue;
1747 if (r>250) continue;
1748 for (Double_t z= -250; z<250; z+=step){
1749 Int_t roc=(z>0)?0:18;
1750 xyz[0]=x;
1751 xyz[1]=y;
1752 xyz[2]=z;
1753 Double_t phi = TMath::ATan2(y,x);
1754 DistortPoint(xyz,roc);
1755 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1756 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
1757 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1758 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1759 Double_t dx = xyz[0]-x;
1760 Double_t dy = xyz[1]-y;
1761 Double_t dz = xyz[2]-z;
1762 Double_t dr=r1-r;
1763 Double_t drphi=(phi1-phi)*r;
1764 (*pcstream)<<"distortion"<<
1765 "x="<<x<< // original position
1766 "y="<<y<<
1767 "z="<<z<<
1768 "r="<<r<<
1769 "phi="<<phi<<
1770 "x1="<<xyz[0]<< // distorted position
1771 "y1="<<xyz[1]<<
1772 "z1="<<xyz[2]<<
1773 "r1="<<r1<<
1774 "phi1="<<phi1<<
1775 //
1776 "dx="<<dx<< // delta position
1777 "dy="<<dy<<
1778 "dz="<<dz<<
1779 "dr="<<dr<<
1780 "drphi="<<drphi<<
1781 "\n";
1782 }
1783 }
1784 }
1785 delete pcstream;
1786 TFile f(Form("correction%s.root",GetName()));
1787 TTree * tree = (TTree*)f.Get("distortion");
1788 TTree * tree2= tree->CopyTree("1");
1789 tree2->SetName(Form("dist%s",GetName()));
1790 tree2->SetDirectory(0);
1791 delete tree;
1792 return tree2;
1793}
1794
1795
1796
be67055b 1797
46e89793 1798void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
be67055b 1799 //
1800 // Make a fit tree:
1801 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1802 // calculates partial distortions
1803 // Partial distortion is stored in the resulting tree
1804 // Output is storred in the file distortion_<dettype>_<partype>.root
1805 // Partial distortion is stored with the name given by correction name
1806 //
1807 //
1808 // Parameters of function:
1809 // input - input tree
cfe2c39a 1810 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing
be67055b 1811 // ppype - parameter type
1812 // corrArray - array with partial corrections
1813 // step - skipe entries - if 1 all entries processed - it is slow
1814 // debug 0 if debug on also space points dumped - it is slow
c9cbd2f2 1815
b322e06a 1816 const Double_t kMaxSnp = 0.85;
cfe2c39a 1817 const Double_t kcutSnp=0.25;
1818 const Double_t kcutTheta=1.;
1819 const Double_t kRadiusTPC=85;
1820 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1821 //
b322e06a 1822 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1823 // const Double_t kB2C=-0.299792458e-3;
46e89793 1824 const Int_t kMinEntries=20;
cfe2c39a 1825 Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
46e89793 1826 Float_t refX;
1827 Int_t run;
1828 tinput->SetBranchAddress("run",&run);
be67055b 1829 tinput->SetBranchAddress("theta",&theta);
1830 tinput->SetBranchAddress("phi", &phi);
1831 tinput->SetBranchAddress("snp",&snp);
1832 tinput->SetBranchAddress("mean",&mean);
1833 tinput->SetBranchAddress("rms",&rms);
1834 tinput->SetBranchAddress("entries",&entries);
cfe2c39a 1835 tinput->SetBranchAddress("sector",&sector);
1836 tinput->SetBranchAddress("dsec",&dsec);
1837 tinput->SetBranchAddress("refX",&refX);
46e89793 1838 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
be67055b 1839 //
1840 Int_t nentries=tinput->GetEntries();
1841 Int_t ncorr=corrArray->GetEntries();
7f4cb119 1842 Double_t corrections[100]={0}; //
be67055b 1843 Double_t tPar[5];
1844 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
be67055b 1845 Int_t dir=0;
cfe2c39a 1846 if (dtype==5 || dtype==6) dtype=4;
1847 if (dtype==0) { dir=-1;}
1848 if (dtype==1) { dir=1;}
1849 if (dtype==2) { dir=-1;}
1850 if (dtype==3) { dir=1;}
1851 if (dtype==4) { dir=-1;}
be67055b 1852 //
46e89793 1853 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
be67055b 1854 tinput->GetEntry(ientry);
7f4cb119 1855 if (TMath::Abs(snp)>kMaxSnp) continue;
be67055b 1856 tPar[0]=0;
1857 tPar[1]=theta*refX;
cfe2c39a 1858 if (dtype==2) tPar[1]=theta*kRadiusTPC;
be67055b 1859 tPar[2]=snp;
1860 tPar[3]=theta;
4486a91f 1861 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
cfe2c39a 1862 if (dtype==4){
1863 // tracks crossing CE
1864 tPar[1]=0; // track at the CE
1865 //if (TMath::Abs(theta) <0.05) continue; // deep cross
1866 }
1867
1868 if (TMath::Abs(snp) >kcutSnp) continue;
1869 if (TMath::Abs(theta) >kcutTheta) continue;
1870 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
8b63d99c 1871 Double_t bz=AliTrackerBase::GetBz();
cfe2c39a 1872 if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks
1873 if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1874
1875 if (dtype==2 && TMath::Abs(bz)>0.1 ) {
1876 tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1877 // snp at the TPC inner radius in case the vertex match used
1878 }
1879 }
1880 //
4486a91f 1881 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
7f4cb119 1882 AliExternalTrackParam track(refX,phi,tPar,cov);
1883 Double_t xyz[3];
1884 track.GetXYZ(xyz);
1885 Int_t id=0;
46e89793 1886 Double_t pt=1./tPar[4];
7f4cb119 1887 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
cfe2c39a 1888 //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
46e89793 1889 Double_t refXD=refX;
be67055b 1890 (*pcstream)<<"fit"<<
46e89793 1891 "run="<<run<< // run number
8b63d99c 1892 "bz="<<bz<< // magnetic filed used
be67055b 1893 "dtype="<<dtype<< // detector match type
1894 "ptype="<<ptype<< // parameter type
1895 "theta="<<theta<< // theta
1896 "phi="<<phi<< // phi
1897 "snp="<<snp<< // snp
1898 "mean="<<mean<< // mean dist value
1899 "rms="<<rms<< // rms
cfe2c39a 1900 "sector="<<sector<<
1901 "dsec="<<dsec<<
46e89793 1902 "refX="<<refXD<< // referece X as double
7f4cb119 1903 "gx="<<xyz[0]<< // global position at reference
1904 "gy="<<xyz[1]<< // global position at reference
1905 "gz="<<xyz[2]<< // global position at reference
1906 "dRrec="<<dRrec<< // delta Radius in reconstruction
46e89793 1907 "pt="<<pt<< // pt
7f4cb119 1908 "id="<<id<< // track id
be67055b 1909 "entries="<<entries;// number of entries in bin
1910 //
cfe2c39a 1911 Bool_t isOK=kTRUE;
46e89793 1912 if (entries<kMinEntries) isOK=kFALSE;
1913 //
cfe2c39a 1914 if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
be67055b 1915 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1916 corrections[icorr]=0;
1917 if (entries>kMinEntries){
1918 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1919 AliExternalTrackParam *trackOut = 0;
1920 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1921 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
cfe2c39a 1922 if (dtype==0) {dir= -1;}
1923 if (dtype==1) {dir= 1;}
1924 if (dtype==2) {dir= -1;}
1925 if (dtype==3) {dir= 1;}
b1f0a2a5 1926 //
7f4cb119 1927 if (trackOut){
cfe2c39a 1928 if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1929 if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
1930 if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1931 // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1932 //
7f4cb119 1933 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1934 delete trackOut;
1935 }else{
1936 corrections[icorr]=0;
cfe2c39a 1937 isOK=kFALSE;
7f4cb119 1938 }
cfe2c39a 1939 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
be67055b 1940 }
1941 (*pcstream)<<"fit"<<
46e89793 1942 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
be67055b 1943 }
cfe2c39a 1944
1945 if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1946 //
1947 // special case of the TPC tracks crossing the CE
1948 //
1949 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1950 corrections[icorr]=0;
1951 if (entries>kMinEntries){
46e89793 1952 AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
1953 AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
cfe2c39a 1954 AliExternalTrackParam *trackOut0 = 0;
1955 AliExternalTrackParam *trackOut1 = 0;
1956 //
1957 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1958 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1959 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1960 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1961 //
1962 if (trackOut0 && trackOut1){
1963 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1964 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1965 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1966 if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1967 //
1968 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1969 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1970 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1971 if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
1972 if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1973 //
1974 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1975 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
46e89793 1976 if (isOK)
1977 if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
1978 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
1979 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
1980 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
1981 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
1982 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
1983 ){
1984 isOK=kFALSE;
1985 }
cfe2c39a 1986 delete trackOut0;
46e89793 1987 delete trackOut1;
cfe2c39a 1988 }else{
1989 corrections[icorr]=0;
1990 isOK=kFALSE;
1991 }
1992 //
1993 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
1994 }
cfe2c39a 1995 (*pcstream)<<"fit"<<
46e89793 1996 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
cfe2c39a 1997 }
1998 //
1999 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
be67055b 2000 }
cfe2c39a 2001
2002
be67055b 2003 delete pcstream;
2004}
2005
2006
2007
46e89793 2008void AliTPCCorrection::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
2009 //
2010 // Make a fit tree:
2011 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
2012 // calculates partial distortions
2013 // Partial distortion is stored in the resulting tree
2014 // Output is storred in the file distortion_<dettype>_<partype>.root
2015 // Partial distortion is stored with the name given by correction name
2016 //
2017 //
2018 // Parameters of function:
2019 // input - input tree
2020 // dtype - distortion type 10 - IROC-OROC
2021 // ppype - parameter type
2022 // corrArray - array with partial corrections
2023 // step - skipe entries - if 1 all entries processed - it is slow
2024 // debug 0 if debug on also space points dumped - it is slow
2025
2026 const Double_t kMaxSnp = 0.8;
2027 const Int_t kMinEntries=200;
2028 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
2029 //
2030 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2031 // const Double_t kB2C=-0.299792458e-3;
2032 Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
2033 Int_t isec1, isec0;
2034 Double_t refXD;
2035 Float_t refX;
2036 Int_t run;
2037 tinput->SetBranchAddress("run",&run);
2038 tinput->SetBranchAddress("theta",&theta);
2039 tinput->SetBranchAddress("phi", &phi);
2040 tinput->SetBranchAddress("snp",&snp);
2041 tinput->SetBranchAddress("mean",&mean);
2042 tinput->SetBranchAddress("rms",&rms);
2043 tinput->SetBranchAddress("entries",&entries);
2044 tinput->SetBranchAddress("sector",&sector);
2045 tinput->SetBranchAddress("dsec",&dsec);
2046 tinput->SetBranchAddress("refX",&refXD);
2047 tinput->SetBranchAddress("z",&globalZ);
2048 tinput->SetBranchAddress("isec0",&isec0);
2049 tinput->SetBranchAddress("isec1",&isec1);
2050 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
2051 //
2052 Int_t nentries=tinput->GetEntries();
2053 Int_t ncorr=corrArray->GetEntries();
2054 Double_t corrections[100]={0}; //
2055 Double_t tPar[5];
2056 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
2057 Int_t dir=0;
2058 //
2059 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
2060 tinput->GetEntry(ientry);
2061 refX=refXD;
2062 Int_t id=-1;
2063 if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
2064 if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
2065 if (dtype==10 && id==-1) continue;
2066 //
2067 dir=-1;
2068 tPar[0]=0;
2069 tPar[1]=globalZ;
2070 tPar[2]=snp;
2071 tPar[3]=theta;
2072 tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
2073 Double_t pt=1./tPar[4];
2074 //
2075 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
2076 Double_t bz=AliTrackerBase::GetBz();
2077 AliExternalTrackParam track(refX,phi,tPar,cov);
2078 Double_t xyz[3],xyzIn[3],xyzOut[3];
2079 track.GetXYZ(xyz);
2080 track.GetXYZAt(85,bz,xyzIn);
2081 track.GetXYZAt(245,bz,xyzOut);
2082 Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
2083 Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
2084 Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
2085 Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
2086 Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
2087 Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
2088 //
2089 Bool_t isOK=kTRUE;
2090 if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
2091 if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
2092 if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
2093 // Do downscale
2094 if (TMath::Abs(theta)>1) isOK=kFALSE;
2095 //
2096 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
2097 //
2098 (*pcstream)<<"fit"<<
2099 "run="<<run<< //run
2100 "bz="<<bz<< // magnetic filed used
2101 "dtype="<<dtype<< // detector match type
2102 "ptype="<<ptype<< // parameter type
2103 "theta="<<theta<< // theta
2104 "phi="<<phi<< // phi
2105 "snp="<<snp<< // snp
2106 "mean="<<mean<< // mean dist value
2107 "rms="<<rms<< // rms
2108 "sector="<<sector<<
2109 "dsec="<<dsec<<
2110 "refX="<<refXD<< // referece X
2111 "gx="<<xyz[0]<< // global position at reference
2112 "gy="<<xyz[1]<< // global position at reference
2113 "gz="<<xyz[2]<< // global position at reference
2114 "dRrec="<<dRrec<< // delta Radius in reconstruction
2115 "pt="<<pt<< //pt
2116 "id="<<id<< // track id
2117 "entries="<<entries;// number of entries in bin
2118 //
2119 AliExternalTrackParam *trackOut0 = 0;
2120 AliExternalTrackParam *trackOut1 = 0;
2121 AliExternalTrackParam *ptrackIn0 = 0;
2122 AliExternalTrackParam *ptrackIn1 = 0;
2123
2124 for (Int_t icorr=0; icorr<ncorr; icorr++) {
2125 //
2126 // special case of the TPC tracks crossing the CE
2127 //
2128 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2129 corrections[icorr]=0;
2130 if (entries>kMinEntries &&isOK){
2131 AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
2132 AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
2133 ptrackIn1=&trackIn0;
2134 ptrackIn0=&trackIn1;
2135 //
2136 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
2137 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
2138 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
2139 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
2140 //
2141 if (trackOut0 && trackOut1){
2142 //
2143 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
2144 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2145 // rotate all tracks to the same frame
2146 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2147 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2148 if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2149 //
2150 if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2151 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2152 if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2153 //
2154 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
2155 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
2156 (*pcstream)<<"fitDebug"<< // just to debug the correction
2157 "mean="<<mean<<
2158 "pIn0.="<<ptrackIn0<<
2159 "pIn1.="<<ptrackIn1<<
2160 "pOut0.="<<trackOut0<<
2161 "pOut1.="<<trackOut1<<
2162 "refX="<<refXD<<
2163 "\n";
2164 delete trackOut0;
2165 delete trackOut1;
2166 }else{
2167 corrections[icorr]=0;
2168 isOK=kFALSE;
2169 }
2170 }
2171 (*pcstream)<<"fit"<<
2172 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
2173 }
2174 //
2175 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2176 }
2177 delete pcstream;
2178}
2179
2180
2181
2182void AliTPCCorrection::MakeLaserDistortionTreeOld(TTree* tree, TObjArray *corrArray, Int_t itype){
7f4cb119 2183 //
2184 // Make a laser fit tree for global minimization
2185 //
2186 const Double_t cutErrY=0.1;
2187 const Double_t cutErrZ=0.1;
2188 const Double_t kEpsilon=0.00000001;
46e89793 2189 const Double_t kMaxDist=1.; // max distance - space correction
2190 const Double_t kMaxRMS=0.05; // max distance -between point and local mean
7f4cb119 2191 TVectorD *vecdY=0;
2192 TVectorD *vecdZ=0;
2193 TVectorD *veceY=0;
2194 TVectorD *veceZ=0;
2195 AliTPCLaserTrack *ltr=0;
2196 AliTPCLaserTrack::LoadTracks();
2197 tree->SetBranchAddress("dY.",&vecdY);
2198 tree->SetBranchAddress("dZ.",&vecdZ);
2199 tree->SetBranchAddress("eY.",&veceY);
2200 tree->SetBranchAddress("eZ.",&veceZ);
2201 tree->SetBranchAddress("LTr.",&ltr);
2202 Int_t entries= tree->GetEntries();
cfe2c39a 2203 TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
7f4cb119 2204 Double_t bz=AliTrackerBase::GetBz();
2205 //
2206
2207 for (Int_t ientry=0; ientry<entries; ientry++){
2208 tree->GetEntry(ientry);
2209 if (!ltr->GetVecGX()){
2210 ltr->UpdatePoints();
2211 }
2212 TVectorD * delta= (itype==0)? vecdY:vecdZ;
2213 TVectorD * err= (itype==0)? veceY:veceZ;
46e89793 2214 TLinearFitter fitter(2,"pol1");
2215 for (Int_t iter=0; iter<2; iter++){
2216 Double_t kfit0=0, kfit1=0;
2217 Int_t npoints=fitter.GetNpoints();
2218 if (npoints>80){
2219 fitter.Eval();
2220 kfit0=fitter.GetParameter(0);
2221 kfit1=fitter.GetParameter(1);
2222 }
2223 for (Int_t irow=0; irow<159; irow++){
2224 Bool_t isOK=kTRUE;
2225 Int_t isOKF=0;
2226 Int_t nentries = 1000;
2227 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
2228 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
2229 Int_t dtype=5;
2230 Double_t array[10];
2231 Int_t first3=TMath::Max(irow-3,0);
2232 Int_t last3 =TMath::Min(irow+3,159);
2233 Int_t counter=0;
2234 if ((*ltr->GetVecSec())[irow]>=0 && err) {
2235 for (Int_t jrow=first3; jrow<=last3; jrow++){
2236 if ((*ltr->GetVecSec())[irow]!= (*ltr->GetVecSec())[jrow]) continue;
2237 if ((*err)[jrow]<kEpsilon) continue;
2238 array[counter]=(*delta)[jrow];
2239 counter++;
2240 }
2241 }
2242 Double_t rms3 = 0;
2243 Double_t mean3 = 0;
2244 if (counter>2){
2245 rms3 = TMath::RMS(counter,array);
2246 mean3 = TMath::Mean(counter,array);
2247 }else{
2248 isOK=kFALSE;
2249 }
2250 Double_t phi =(*ltr->GetVecPhi())[irow];
2251 Double_t theta =ltr->GetTgl();
2252 Double_t mean=delta->GetMatrixArray()[irow];
2253 Double_t gx=0,gy=0,gz=0;
2254 Double_t snp = (*ltr->GetVecP2())[irow];
2255 Double_t dRrec=0;
2256 // Double_t rms = err->GetMatrixArray()[irow];
cfe2c39a 2257 //
46e89793 2258 gx = (*ltr->GetVecGX())[irow];
2259 gy = (*ltr->GetVecGY())[irow];
2260 gz = (*ltr->GetVecGZ())[irow];
2261 //
2262 // get delta R used in reconstruction
2263 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2264 AliTPCCorrection * correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
2265 // const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
2266 //Double_t xyz0[3]={gx,gy,gz};
2267 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
2268 Double_t fphi = TMath::ATan2(gy,gx);
2269 Double_t fsector = 9.*fphi/TMath::Pi();
2270 if (fsector<0) fsector+=18;
2271 Double_t dsec = fsector-Int_t(fsector)-0.5;
2272 Double_t refX=0;
2273 Int_t id= ltr->GetId();
2274 Double_t pt=0;
2275 //
2276 if (1 && oldR>1) {
2277 Float_t xyz1[3]={gx,gy,gz};
2278 Int_t sector=(gz>0)?0:18;
2279 correction->CorrectPoint(xyz1, sector);
2280 refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
2281 dRrec=oldR-refX;
2282 }
2283 if (TMath::Abs(rms3)>kMaxRMS) isOK=kFALSE;
2284 if (TMath::Abs(mean-mean3)>kMaxRMS) isOK=kFALSE;
2285 if (counter<4) isOK=kFALSE;
2286 if (npoints<90) isOK=kFALSE;
2287 if (isOK){
2288 fitter.AddPoint(&refX,mean);
7f4cb119 2289 }
46e89793 2290 Double_t deltaF=kfit0+kfit1*refX;
2291 if (iter==1){
2292 (*pcstream)<<"fitFull"<< // dumpe also intermediate results
2293 "bz="<<bz<< // magnetic filed used
2294 "dtype="<<dtype<< // detector match type
2295 "ptype="<<itype<< // parameter type
2296 "theta="<<theta<< // theta
2297 "phi="<<phi<< // phi
2298 "snp="<<snp<< // snp
2299 "mean="<<mean3<< // mean dist value
2300 "rms="<<rms3<< // rms
2301 "deltaF="<<deltaF<<
2302 "npoints="<<npoints<< //number of points
2303 "mean3="<<mean3<< // mean dist value
2304 "rms3="<<rms3<< // rms
2305 "counter="<<counter<<
2306 "sector="<<fsector<<
2307 "dsec="<<dsec<<
2308 //
2309 "refX="<<refX<< // reference radius
2310 "gx="<<gx<< // global position
2311 "gy="<<gy<< // global position
2312 "gz="<<gz<< // global position
2313 "dRrec="<<dRrec<< // delta Radius in reconstruction
2314 "id="<<id<< //bundle
2315 "entries="<<nentries<<// number of entries in bin
2316 "\n";
2317 }
2318 if (iter==1) (*pcstream)<<"fit"<< // dump valus for fit
2319 "bz="<<bz<< // magnetic filed used
2320 "dtype="<<dtype<< // detector match type
2321 "ptype="<<itype<< // parameter type
2322 "theta="<<theta<< // theta
2323 "phi="<<phi<< // phi
2324 "snp="<<snp<< // snp
2325 "mean="<<mean3<< // mean dist value
2326 "rms="<<rms3<< // rms
2327 "sector="<<fsector<<
2328 "dsec="<<dsec<<
2329 //
2330 "refX="<<refX<< // reference radius
2331 "gx="<<gx<< // global position
2332 "gy="<<gy<< // global position
2333 "gz="<<gz<< // global position
2334 "dRrec="<<dRrec<< // delta Radius in reconstruction
2335 "pt="<<pt<< //pt
2336 "id="<<id<< //bundle
2337 "entries="<<nentries;// number of entries in bin
2338 //
2339 //
2340 Double_t ky = TMath::Tan(TMath::ASin(snp));
2341 Int_t ncorr = corrArray->GetEntries();
2342 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
2343 Double_t phi0 = TMath::ATan2(gy,gx);
2344 Double_t distortions[1000]={0};
2345 Double_t distortionsR[1000]={0};
2346 if (iter==1){
2347 for (Int_t icorr=0; icorr<ncorr; icorr++) {
2348 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2349 Float_t distPoint[3]={gx,gy,gz};
2350 Int_t sector= (gz>0)? 0:18;
2351 if (r0>80){
2352 corr->DistortPoint(distPoint, sector);
2353 }
2354 // Double_t value=distPoint[2]-gz;
2355 if (itype==0 && r0>1){
2356 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2357 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
2358 Double_t drphi= r0*(phi1-phi0);
2359 Double_t dr = r1-r0;
2360 distortions[icorr] = drphi-ky*dr;
2361 distortionsR[icorr] = dr;
2362 }
2363 if (TMath::Abs(distortions[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE; }
2364 if (TMath::Abs(distortionsR[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE;}
2365 (*pcstream)<<"fit"<<
2366 Form("%s=",corr->GetName())<<distortions[icorr]; // dump correction value
2367 }
2368 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
7f4cb119 2369 }
7f4cb119 2370 }
7f4cb119 2371 }
2372 }
2373 delete pcstream;
2374}
2375
2376
be67055b 2377
97d17739 2378void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ){
cfe2c39a 2379 //
2380 // make a distortion map out ou fthe residual histogram
2381 // Results are written to the debug streamer - pcstream
2382 // Parameters:
2383 // his0 - input (4D) residual histogram
2384 // pcstream - file to write the tree
2385 // run - run number
2386 // refX - track matching reference X
2387 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
2388 // THnSparse axes:
2389 // OBJ: TAxis #Delta #Delta
2390 // OBJ: TAxis tanTheta tan(#Theta)
2391 // OBJ: TAxis phi #phi
2392 // OBJ: TAxis snp snp
2393
2394 // marian.ivanov@cern.ch
2395 const Int_t kMinEntries=10;
2396 Double_t bz=AliTrackerBase::GetBz();
2397 Int_t idim[4]={0,1,2,3};
2398 //
2399 //
2400 //
2401 Int_t nbins3=his0->GetAxis(3)->GetNbins();
2402 Int_t first3=his0->GetAxis(3)->GetFirst();
2403 Int_t last3 =his0->GetAxis(3)->GetLast();
2404 //
2405 for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle
97d17739 2406 his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
cfe2c39a 2407 Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
2408 THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
2409 //
2410 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
2411 Int_t first2 = his3->GetAxis(2)->GetFirst();
2412 Int_t last2 = his3->GetAxis(2)->GetLast();
2413 //
2414 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
97d17739 2415 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
cfe2c39a 2416 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
2417 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
2418 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
2419 Int_t first1 = his2->GetAxis(1)->GetFirst();
2420 Int_t last1 = his2->GetAxis(1)->GetLast();
2421 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
2422 //
2423 Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
2424 his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
2425 if (TMath::Abs(x1)<0.1){
2426 if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
2427 if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
2428 }
2429 if (TMath::Abs(x1)<0.06){
2430 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
2431 }
2432 TH1 * hisDelta = his2->Projection(0);
2433 //
2434 Double_t entries = hisDelta->GetEntries();
2435 Double_t mean=0, rms=0;
2436 if (entries>kMinEntries){
2437 mean = hisDelta->GetMean();
2438 rms = hisDelta->GetRMS();
2439 }
2440 Double_t sector = 9.*x2/TMath::Pi();
2441 if (sector<0) sector+=18;
2442 Double_t dsec = sector-Int_t(sector)-0.5;
2443 Double_t z=refX*x1;
2444 (*pcstream)<<hname<<
2445 "run="<<run<<
2446 "bz="<<bz<<
2447 "theta="<<x1<<
2448 "phi="<<x2<<
2449 "z="<<z<< // dummy z
2450 "snp="<<x3<<
2451 "entries="<<entries<<
2452 "mean="<<mean<<
2453 "rms="<<rms<<
2454 "refX="<<refX<< // track matching refernce plane
2455 "type="<<type<< //
2456 "sector="<<sector<<
2457 "dsec="<<dsec<<
2458 "\n";
2459 delete hisDelta;
02cd5ade 2460 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
cfe2c39a 2461 }
2462 delete his2;
2463 }
2464 delete his3;
2465 }
2466}
2467
2468
2469
2470
2471void AliTPCCorrection::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
8b63d99c 2472 //
2473 // make a distortion map out ou fthe residual histogram
2474 // Results are written to the debug streamer - pcstream
2475 // Parameters:
2476 // his0 - input (4D) residual histogram
2477 // pcstream - file to write the tree
2478 // run - run number
cfe2c39a 2479 // refX - track matching reference X
2480 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
8b63d99c 2481 // marian.ivanov@cern.ch
cfe2c39a 2482 //
2483 // Histo axeses
2484 // Collection name='TObjArray', class='TObjArray', size=16
2485 // 0. OBJ: TAxis #Delta #Delta
2486 // 1. OBJ: TAxis N_{cl} N_{cl}
2487 // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm)
2488 // 3. OBJ: TAxis z (cm) z (cm)
2489 // 4. OBJ: TAxis sin(#phi) sin(#phi)
2490 // 5. OBJ: TAxis tan(#theta) tan(#theta)
2491 // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV)
2492 // 7. OBJ: TAxis pt (GeV) pt (GeV)
2493 // 8. OBJ: TAxis alpha alpha
2494 const Int_t kMinEntries=10;
2495 //
2496 // 1. make default selections
2497 //
2498 TH1 * hisDelta=0;
2499 Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z
2500 hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks
2501 hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe
2502 hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
2503 hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
2504 hisDelta= hisInput->Projection(0);
2505 hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
2506 delete hisDelta;
2507 THnSparse *his0= hisInput->Projection(4,idim0);
2508 //
2509 // 2. Get mean in diferent bins
2510 //
8b63d99c 2511 Int_t nbins1=his0->GetAxis(1)->GetNbins();
2512 Int_t first1=his0->GetAxis(1)->GetFirst();
2513 Int_t last1 =his0->GetAxis(1)->GetLast();
2514 //
2515 Double_t bz=AliTrackerBase::GetBz();
cfe2c39a 2516 Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
2517 //
2518 for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
2519 //
2520 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
2521 his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
8b63d99c 2522 //
8b63d99c 2523 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
2524 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
2525 Int_t first3 = his1->GetAxis(3)->GetFirst();
2526 Int_t last3 = his1->GetAxis(3)->GetLast();
2527 //
cfe2c39a 2528 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex"
8b63d99c 2529 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
2530 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
2531 if (ibin3<first3) {
2532 his1->GetAxis(3)->SetRangeUser(-1,1);
2533 x3=0;
2534 }
2535 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
2536 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
2537 Int_t first2 = his3->GetAxis(2)->GetFirst();
2538 Int_t last2 = his3->GetAxis(2)->GetLast();
2539 //
cfe2c39a 2540 for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
8b63d99c 2541 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
2542 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
cfe2c39a 2543 hisDelta = his3->Projection(0);
8b63d99c 2544 //
2545 Double_t entries = hisDelta->GetEntries();
2546 Double_t mean=0, rms=0;
2547 if (entries>kMinEntries){
2548 mean = hisDelta->GetMean();
2549 rms = hisDelta->GetRMS();
2550 }
cfe2c39a 2551 Double_t sector = 9.*x2/TMath::Pi();
2552 if (sector<0) sector+=18;
2553 Double_t dsec = sector-Int_t(sector)-0.5;
2554 Double_t snp=0; // dummy snp - equal 0
8b63d99c 2555 (*pcstream)<<hname<<
2556 "run="<<run<<
cfe2c39a 2557 "bz="<<bz<< // magnetic field
2558 "theta="<<x1<< // theta
2559 "phi="<<x2<< // phi (alpha)
2560 "z="<<x3<< // z at "vertex"
2561 "snp="<<snp<< // dummy snp
2562 "entries="<<entries<< // entries in bin
2563 "mean="<<mean<< // mean
8b63d99c 2564 "rms="<<rms<<
cfe2c39a 2565 "refX="<<refX<< // track matching refernce plane
2566 "type="<<type<< // parameter type
2567 "sector="<<sector<< // sector
2568 "dsec="<<dsec<< // dummy delta sector
8b63d99c 2569 "\n";
2570 delete hisDelta;
2571 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
2572 }
2573 delete his3;
2574 }
2575 delete his1;
2576 }
cfe2c39a 2577 delete his0;
2578}
2579
2580
2581
2582void AliTPCCorrection::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type){
2583 //
2584 // make a distortion map out of the residual histogram
2585 // Results are written to the debug streamer - pcstream
2586 // Parameters:
2587 // his0 - input (4D) residual histogram
2588 // pcstream - file to write the tree
2589 // run - run number
2590 // type - 0- y 1-z,2 -snp, 3-theta
2591 // marian.ivanov@cern.ch
2592
2593 //Collection name='TObjArray', class='TObjArray', size=16
2594 //0 OBJ: TAxis delta delta
2595 //1 OBJ: TAxis phi phi
2596 //2 OBJ: TAxis localX localX
2597 //3 OBJ: TAxis kY kY
2598 //4 OBJ: TAxis kZ kZ
2599 //5 OBJ: TAxis is1 is1
2600 //6 OBJ: TAxis is0 is0
2601 //7. OBJ: TAxis z z
2602 //8. OBJ: TAxis IsPrimary IsPrimary
2603
2604 const Int_t kMinEntries=10;
2605 THnSparse * hisSector0=0;
2606 TH1 * htemp=0; // histogram to calculate mean value of parameter
2607 Double_t bz=AliTrackerBase::GetBz();
2608
2609 //
2610 // Loop over pair of sector:
2611 // isPrim - 8 ==> 8
2612 // isec0 - 6 ==> 7
2613 // isec1 - 5 ==> 6
2614 // refX - 2 ==> 5
2615 //
2616 // phi - 1 ==> 4
2617 // z - 7 ==> 3
2618 // snp - 3 ==> 2
2619 // theta- 4 ==> 1
2620 // 0 ==> 0;
2621 for (Int_t isec0=0; isec0<72; isec0++){
2622 Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
2623 //
2624 //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ?
2625 hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
2626 hisSector0=hisInput->Projection(7,index0);
2627 //
2628 //
2629 for (Int_t isec1=isec0+1; isec1<72; isec1++){
2630 //if (isec1!=isec0+36) continue;
2631 if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
2632 printf("Sectors %d\t%d\n",isec1,isec0);
2633 hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
2634 TH1 * hisX=hisSector0->Projection(5);
2635 Double_t refX= hisX->GetMean();
2636 delete hisX;
2637 TH1 *hisDelta=hisSector0->Projection(0);
2638 Double_t dmean = hisDelta->GetMean();
2639 Double_t drms = hisDelta->GetRMS();
2640 hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
2641 delete hisDelta;
2642 //
2643 // 1. make default selections
2644 //
2645 Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
2646 THnSparse *hisSector1= hisSector0->Projection(5,idim0);
2647 //
2648 // 2. Get mean in diferent bins
2649 //
2650 Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4}
2651 //
2652 // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
2653 Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
2654 Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
2655 //
2656 for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=1){ //axis 4 - phi
2657 //
2658 // Phi loop
2659 //
2660 Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
2661 Double_t psec = (9*xPhi/TMath::Pi());
2662 if (psec<0) psec+=18;
2663 Bool_t isOK0=kFALSE;
2664 Bool_t isOK1=kFALSE;
2665 if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE;
2666 if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE;
2667 if (!isOK0) continue;
2668 if (!isOK1) continue;
2669 //
2670 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
2671 if (isec1!=isec0+36) {
2672 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
2673 }
2674 //
2675 htemp = hisSector1->Projection(4);
2676 xPhi=htemp->GetMean();
2677 delete htemp;
2678 THnSparse * hisPhi = hisSector1->Projection(4,idim);
2679 //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins();
2680 Int_t firstZ = hisPhi->GetAxis(3)->GetFirst();
2681 Int_t lastZ = hisPhi->GetAxis(3)->GetLast();
2682 //
2683 for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=1){ // axis 3 - z
2684 //
2685 // Z loop
2686 //
2687 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
2688 if (isec1!=isec0+36) {
2689 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
2690 }
2691 htemp = hisPhi->Projection(3);
2692 Double_t xZ= htemp->GetMean();
2693 delete htemp;
2694 THnSparse * hisZ= hisPhi->Projection(3,idim);
2695 //projected histogram according selection 3 -z
2696 //
2697 //
2698 //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins();
2699 Int_t firstSnp = hisZ->GetAxis(2)->GetFirst();
2700 Int_t lastSnp = hisZ->GetAxis(2)->GetLast();
2701 for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp
2702 //
2703 // Snp loop
2704 //
2705 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
2706 if (isec1!=isec0+36) {
2707 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
2708 }
2709 htemp = hisZ->Projection(2);
2710 Double_t xSnp= htemp->GetMean();
2711 delete htemp;
2712 THnSparse * hisSnp= hisZ->Projection(2,idim);
2713 //projected histogram according selection 2 - snp
2714
2715 //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
2716 Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
2717 Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
2718 //
2719 for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
2720
2721
2722 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
2723 if (isec1!=isec0+36) {
2724 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
2725 }
2726 htemp = hisSnp->Projection(1);
2727 Double_t xTheta=htemp->GetMean();
2728 delete htemp;
2729 hisDelta = hisSnp->Projection(0);
2730 //
2731 Double_t entries = hisDelta->GetEntries();
2732 Double_t mean=0, rms=0;
2733 if (entries>kMinEntries){
2734 mean = hisDelta->GetMean();
2735 rms = hisDelta->GetRMS();
2736 }
2737 Double_t sector = 9.*xPhi/TMath::Pi();
2738 if (sector<0) sector+=18;
2739 Double_t dsec = sector-Int_t(sector)-0.5;
2740 Int_t dtype=1; // TPC alignment type
2741 (*pcstream)<<hname<<
2742 "run="<<run<<
2743 "bz="<<bz<< // magnetic field
2744 "ptype="<<type<< // parameter type
2745 "dtype="<<dtype<< // parameter type
2746 "isec0="<<isec0<< // sector 0
2747 "isec1="<<isec1<< // sector 1
2748 "sector="<<sector<< // sector as float
2749 "dsec="<<dsec<< // delta sector
2750 //
2751 "theta="<<xTheta<< // theta
2752 "phi="<<xPhi<< // phi (alpha)
2753 "z="<<xZ<< // z
2754 "snp="<<xSnp<< // snp
2755 //
2756 "entries="<<entries<< // entries in bin
2757 "mean="<<mean<< // mean
2758 "rms="<<rms<< // rms
2759 "refX="<<refX<< // track matching reference plane
2760 "\n";
2761 delete hisDelta;
2762 printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean);
2763 //
2764 }//ibinTheta
2765 delete hisSnp;
2766 } //ibinSnp
2767 delete hisZ;
2768 }//ibinZ
2769 delete hisPhi;
2770 }//ibinPhi
2771 delete hisSector1;
2772 }//isec1
2773 delete hisSector0;
2774 }//isec0
8b63d99c 2775}
2776
2777
2778
2779
2780
cfe2c39a 2781
2782
ffab0c37 2783void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
2784 //
2785 // Store object in the OCDB
2786 // By default the object is stored in the current directory
2787 // default comment consit of user name and the date
2788 //
2789 TString ocdbStorage="";
2790 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
2791 AliCDBMetaData *metaData= new AliCDBMetaData();
2792 metaData->SetObjectClassName("AliTPCCorrection");
2793 metaData->SetResponsible("Marian Ivanov");
2794 metaData->SetBeamPeriod(1);
2795 metaData->SetAliRootVersion("05-25-01"); //root version
2796 TString userName=gSystem->GetFromPipe("echo $USER");
2797 TString date=gSystem->GetFromPipe("date");
2798
2799 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
2800 if (comment) metaData->SetComment(comment);
2801 AliCDBId* id1=NULL;
2802 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
2803 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
2804 gStorage->Put(this, (*id1), metaData);
2805}
2806
ca58ed4e 2807
7d85e147 2808void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
c9cbd2f2 2809 //
2810 // Fast method to simulate the influence of the given distortion on the vertex reconstruction
2811 //
ca58ed4e 2812
c9cbd2f2 2813 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2814 if (!magF) AliError("Magneticd field - not initialized");
2815 Double_t bz = magF->SolenoidField(); //field in kGauss
9f3b99e2 2816 printf("bz: %f\n",bz);
c9cbd2f2 2817 AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
ca58ed4e 2818
c9cbd2f2 2819 TObjArray aTrk; // Original Track array of Aside
2820 TObjArray daTrk; // Distorted Track array of A side
2821 UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
2822 TObjArray cTrk;
2823 TObjArray dcTrk;
2824 UShort_t *cId = new UShort_t [nTracks];
2825 Int_t id=0;
ca58ed4e 2826 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
7d85e147 2827 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
ca58ed4e 2828 fpt.SetParameters(7.24,0.120);
2829 fpt.SetNpx(10000);
2830 for(Int_t nt=0; nt<nTracks; nt++){
2831 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
7d85e147 2832 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
c9cbd2f2 2833 Double_t pt = fpt.GetRandom(); // momentum for f1
2834 // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
ca58ed4e 2835 Short_t sign=1;
2836 if(gRandom->Rndm() < 0.5){
2837 sign =1;
2838 }else{
2839 sign=-1;
2840 }
2841
2842 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
2843 Double_t pxyz[3];
2844 pxyz[0]=pt*TMath::Cos(phi);
2845 pxyz[1]=pt*TMath::Sin(phi);
2846 pxyz[2]=pt*TMath::Tan(theta);
2847 Double_t cv[21]={0};
2848 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
2849
2850 Double_t refX=1.;
2851 Int_t dir=-1;
2852 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
2853 if (!td) continue;
2854 if (pcstream) (*pcstream)<<"track"<<
2855 "eta="<<eta<<
2856 "theta="<<theta<<
2857 "tOrig.="<<t<<
2858 "td.="<<td<<
2859 "\n";
7d85e147 2860 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
ca58ed4e 2861 if (td){
c9cbd2f2 2862 daTrk.AddLast(td);
2863 aTrk.AddLast(t);
2864 Int_t nn=aTrk.GetEntriesFast();
2865 aId[nn]=id;
ca58ed4e 2866 }
7d85e147 2867 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
ca58ed4e 2868 if (td){
c9cbd2f2 2869 dcTrk.AddLast(td);
2870 cTrk.AddLast(t);
2871 Int_t nn=cTrk.GetEntriesFast();
2872 cId[nn]=id;
ca58ed4e 2873 }
2874 }
c9cbd2f2 2875 id++;
ca58ed4e 2876 }// end of track loop
2877
2878 vertexer->SetTPCMode();
2879 vertexer->SetConstraintOff();
2880
c9cbd2f2 2881 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
2882 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
2883 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
2884 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
ca58ed4e 2885 if (pcstream) (*pcstream)<<"vertex"<<
2886 "x="<<orgVertex[0]<<
2887 "y="<<orgVertex[1]<<
2888 "z="<<orgVertex[2]<<
2889 "av.="<<&aV<< // distorted vertex A side
2890 "cv.="<<&cV<< // distroted vertex C side
2891 "avO.="<<&avOrg<< // original vertex A side
2892 "cvO.="<<&cvOrg<<
2893 "\n";
c9cbd2f2 2894 delete []aId;
2895 delete []cId;
ca58ed4e 2896}
f1817479 2897
2898void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
2899 //
2900 // make correction available for visualization using
2901 // TFormula, TFX and TTree::Draw
2902 // important in order to check corrections and also compute dervied variables
2903 // e.g correction partial derivatives
2904 //
2905 // NOTE - class is not owner of correction
2906 //
cfe2c39a 2907 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray(10000);
2908 if (position>=fgVisualCorrection->GetEntriesFast())
2909 fgVisualCorrection->Expand((position+10)*2);
f1817479 2910 fgVisualCorrection->AddAt(corr, position);
2911}
2912
2913
2914
287fbdfa 2915Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
f1817479 2916 //
2917 // calculate the correction at given position - check the geffCorr
2918 //
cfe2c39a 2919 // corrType return values
2920 // 0 - delta R
2921 // 1 - delta RPhi
2922 // 2 - delta Z
2923 // 3 - delta RPHI
2924 //
f1817479 2925 if (!fgVisualCorrection) return 0;
2926 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
2927 if (!corr) return 0;
25732bff 2928
f1817479 2929 Double_t phi=sector*TMath::Pi()/9.;
287fbdfa 2930 Double_t gx = r*TMath::Cos(phi);
2931 Double_t gy = r*TMath::Sin(phi);
2932 Double_t gz = r*kZ;
f1817479 2933 Int_t nsector=(gz>0) ? 0:18;
2934 //
2935 //
2936 //
2937 Float_t distPoint[3]={gx,gy,gz};
2938 corr->DistortPoint(distPoint, nsector);
2939 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
2940 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2941 Double_t phi0=TMath::ATan2(gy,gx);
2942 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
2943 if (axisType==0) return r1-r0;
2944 if (axisType==1) return (phi1-phi0)*r0;
2945 if (axisType==2) return distPoint[2]-gz;
cfe2c39a 2946 if (axisType==3) return (TMath::Cos(phi)*(distPoint[0]-gx)+ TMath::Cos(phi)*(distPoint[1]-gy));
f1817479 2947 return phi1-phi0;
2948}
2949
2950Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
2951 //
2952 // return correction at given x,y,z
2953 //
2954 if (!fgVisualCorrection) return 0;
2955 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
2956 if (!corr) return 0;
2957 Double_t phi0= TMath::ATan2(gy,gx);
2958 Int_t nsector=(gz>0) ? 0:18;
2959 Float_t distPoint[3]={gx,gy,gz};
2960 corr->DistortPoint(distPoint, nsector);
2961 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
2962 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2963 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
2964 if (axisType==0) return r1-r0;
2965 if (axisType==1) return (phi1-phi0)*r0;
2966 if (axisType==2) return distPoint[2]-gz;
2967 return phi1-phi0;
2968}
46e89793 2969
fdbbc146 2970Double_t AliTPCCorrection::GetCorrXYZDz(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
2971 //
2972 // return correction at given x,y,z
2973 //
2974 if (!fgVisualCorrection) return 0;
2975 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
2976 if (!corr) return 0;
2977 Double_t phi0= TMath::ATan2(gy,gx);
2978 Int_t nsector=(gz>0) ? 0:18;
2979 Float_t distPoint[3]={gx,gy,gz};
2980 Float_t dxyz[3]={gx,gy,gz};
2981 //
2982 corr->GetCorrectionDz(distPoint, nsector,dxyz,delta);
2983 distPoint[0]-=dxyz[0];
2984 distPoint[1]-=dxyz[1];
2985 distPoint[2]-=dxyz[2];
2986 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
2987 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2988 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
2989 if (axisType==0) return r1-r0;
2990 if (axisType==1) return (phi1-phi0)*r0;
2991 if (axisType==2) return distPoint[2]-gz;
2992 return phi1-phi0;
2993}
2994
2995Double_t AliTPCCorrection::GetCorrXYZIntegrateZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
2996 //
2997 // return correction at given x,y,z
2998 //
2999 if (!fgVisualCorrection) return 0;
3000 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3001 if (!corr) return 0;
3002 Double_t phi0= TMath::ATan2(gy,gx);
3003 Int_t nsector=(gz>0) ? 0:18;
3004 Float_t distPoint[3]={gx,gy,gz};
3005 Float_t dxyz[3]={gx,gy,gz};
3006 //
3007 corr->GetCorrectionIntegralDz(distPoint, nsector,dxyz,delta);
3008 distPoint[0]-=dxyz[0];
3009 distPoint[1]-=dxyz[1];
3010 distPoint[2]-=dxyz[2];
3011 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3012 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3013 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3014 if (axisType==0) return r1-r0;
3015 if (axisType==1) return (phi1-phi0)*r0;
3016 if (axisType==2) return distPoint[2]-gz;
3017 return phi1-phi0;
3018}
3019
46e89793 3020
3021
3022
3023
284418bc 3024void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
46e89793 3025 //
3026 // Make a laser fit tree for global minimization
3027 //
3028 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
3029 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
3030 if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
3031 correction->AddVisualCorrection(correction,0); //register correction
3032
284418bc 3033 // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
3034 //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
46e89793 3035 //
3036 const Double_t cutErrY=0.05;
3037 const Double_t kSigmaCut=4;
3038 // const Double_t cutErrZ=0.03;
3039 const Double_t kEpsilon=0.00000001;
284418bc 3040 // const Double_t kMaxDist=1.; // max distance - space correction
46e89793 3041 TVectorD *vecdY=0;
3042 TVectorD *vecdZ=0;
3043 TVectorD *veceY=0;
3044 TVectorD *veceZ=0;
3045 AliTPCLaserTrack *ltr=0;
3046 AliTPCLaserTrack::LoadTracks();
3047 tree->SetBranchAddress("dY.",&vecdY);
3048 tree->SetBranchAddress("dZ.",&vecdZ);
3049 tree->SetBranchAddress("eY.",&veceY);
3050 tree->SetBranchAddress("eZ.",&veceZ);
3051 tree->SetBranchAddress("LTr.",&ltr);
3052 Int_t entries= tree->GetEntries();
3053 TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
3054 Double_t bz=AliTrackerBase::GetBz();
3055 //
284418bc 3056 // Double_t globalXYZ[3];
3057 //Double_t globalXYZCorr[3];
46e89793 3058 for (Int_t ientry=0; ientry<entries; ientry++){
3059 tree->GetEntry(ientry);
3060 if (!ltr->GetVecGX()){
3061 ltr->UpdatePoints();
3062 }
3063 //
3064 TVectorD fit10(5);
3065 TVectorD fit5(5);
3066 printf("Entry\t%d\n",ientry);
3067 for (Int_t irow0=0; irow0<158; irow0+=1){
3068 //
3069 TLinearFitter fitter10(4,"hyp3");
3070 TLinearFitter fitter5(2,"hyp1");
3071 Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
3072 if (sector<0) continue;
3073 //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
3074
3075 Double_t refX= (*ltr->GetVecLX())[irow0];
3076 Int_t firstRow1 = TMath::Max(irow0-10,0);
3077 Int_t lastRow1 = TMath::Min(irow0+10,158);
3078 Double_t padWidth=(irow0<64)?0.4:0.6;
3079 // make long range fit
3080 for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
3081 if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
3082 if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
3083 if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
3084 Double_t idealX= (*ltr->GetVecLX())[irow1];
3085 Double_t idealY= (*ltr->GetVecLY())[irow1];
284418bc 3086 // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
46e89793 3087 Double_t gx= (*ltr->GetVecGX())[irow1];
3088 Double_t gy= (*ltr->GetVecGY())[irow1];
3089 Double_t gz= (*ltr->GetVecGZ())[irow1];
3090 Double_t measY=(*vecdY)[irow1]+idealY;
3091 Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
3092 // deltaR = R distorted -R ideal
3093 Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
3094 fitter10.AddPoint(xxx,measY,1);
3095 }
3096 Bool_t isOK=kTRUE;
3097 Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
3098 Double_t mean10 =0;// fitter10.GetParameter(0);
3099 Double_t slope10 =0;// fitter10.GetParameter(0);
3100 Double_t cosPart10 = 0;// fitter10.GetParameter(2);
3101 Double_t sinPart10 =0;// fitter10.GetParameter(3);
3102
3103 if (fitter10.GetNpoints()>10){
3104 fitter10.Eval();
3105 rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
3106 mean10 = fitter10.GetParameter(0);
3107 slope10 = fitter10.GetParameter(1);
3108 cosPart10 = fitter10.GetParameter(2);
3109 sinPart10 = fitter10.GetParameter(3);
3110 //
3111 // make short range fit
3112 //
3113 for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
3114 if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
3115 if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
3116 if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
3117 Double_t idealX= (*ltr->GetVecLX())[irow1];
3118 Double_t idealY= (*ltr->GetVecLY())[irow1];
284418bc 3119 // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
46e89793 3120 Double_t gx= (*ltr->GetVecGX())[irow1];
3121 Double_t gy= (*ltr->GetVecGY())[irow1];
3122 Double_t gz= (*ltr->GetVecGZ())[irow1];
3123 Double_t measY=(*vecdY)[irow1]+idealY;
3124 Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
3125 // deltaR = R distorted -R ideal
3126 Double_t expY= mean10+slope10*(idealX+deltaR-refX);
3127 if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
3128 //
3129 Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
3130 Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
3131 fitter5.AddPoint(xxx,measY-corr,1);
3132 }
3133 }else{
3134 isOK=kFALSE;
3135 }
3136 if (fitter5.GetNpoints()<8) isOK=kFALSE;
3137
3138 Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
3139 Double_t offset5 =0;// fitter5.GetParameter(0);
3140 Double_t slope5 =0;// fitter5.GetParameter(0);
3141 if (isOK){
3142 fitter5.Eval();
3143 rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
3144 offset5 = fitter5.GetParameter(0);
3145 slope5 = fitter5.GetParameter(0);
3146 }
3147 //
3148 Double_t dtype=5;
3149 Double_t ptype=0;
3150 Double_t phi =(*ltr->GetVecPhi())[irow0];
3151 Double_t theta =ltr->GetTgl();
3152 Double_t mean=(vecdY)->GetMatrixArray()[irow0];
3153 Double_t gx=0,gy=0,gz=0;
3154 Double_t snp = (*ltr->GetVecP2())[irow0];
3155 Int_t bundle= ltr->GetBundle();
3156 Int_t id= ltr->GetId();
3157 // Double_t rms = err->GetMatrixArray()[irow];
3158 //
3159 gx = (*ltr->GetVecGX())[irow0];
3160 gy = (*ltr->GetVecGY())[irow0];
3161 gz = (*ltr->GetVecGZ())[irow0];
3162 Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
3163 fitter10.GetParameters(fit10);
3164 fitter5.GetParameters(fit5);
3165 Double_t idealY= (*ltr->GetVecLY())[irow0];
3166 Double_t measY=(*vecdY)[irow0]+idealY;
3167 Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
3168 if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
3169 //
3170 (*pcstream)<<"fitFull"<< // dumpe also intermediate results
3171 "bz="<<bz<< // magnetic filed used
3172 "dtype="<<dtype<< // detector match type
3173 "ptype="<<ptype<< // parameter type
3174 "theta="<<theta<< // theta
3175 "phi="<<phi<< // phi
3176 "snp="<<snp<< // snp
3177 "sector="<<sector<<
3178 "bundle="<<bundle<<
3179// // "dsec="<<dsec<<
3180 "refX="<<refX<< // reference radius
3181 "gx="<<gx<< // global position
3182 "gy="<<gy<< // global position
3183 "gz="<<gz<< // global position
3184 "dRrec="<<dRrec<< // delta Radius in reconstruction
3185 "id="<<id<< //bundle
3186 "rms10="<<rms10<<
3187 "rms5="<<rms5<<
3188 "fit10.="<<&fit10<<
3189 "fit5.="<<&fit5<<
3190 "measY="<<measY<<
3191 "mean="<<mean<<
3192 "idealY="<<idealY<<
3193 "corr="<<corr<<
3194 "isOK="<<isOK<<
3195 "\n";
3196 }
3197 }
3198 delete pcstream;
3199}