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