Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /// \class AliTPCTransform
17 : /// \brief Implementation of the TPC transformation class
18 : ///
19 : /// Class for tranformation of the coordinate frame
20 : /// Transformation
21 : /// local coordinate frame (sector, padrow, pad, timebine) ==>
22 : /// rotated global (tracking) cooridnate frame (sector, lx,ly,lz)
23 : ///
24 : /// Unisochronity - (substract time0 - pad by pad)
25 : /// Drift velocity - Currently common drift velocity - functionality of AliTPCParam
26 : /// ExB effect -
27 : ///
28 : /// Time of flight correction -
29 : /// - Depends on the vertex position
30 : /// - by default
31 : ///
32 : /// Usage:
33 : ///
34 : /// ~~~{.cxx}
35 : /// AliTPCclusterer::AddCluster
36 : /// AliTPCtracker::Transform
37 : /// ~~~
38 : ///
39 : /// To test it:
40 : ///
41 : /// ~~~{.cxx}
42 : /// cdb=AliCDBManager::Instance()
43 : /// cdb->SetDefaultStorage("local:///u/mmager/mycalib1")
44 : /// c=AliTPCcalibDB::Instance()
45 : /// c->SetRun(0)
46 : /// Double_t x[]={1.0,2.0,3.0}
47 : /// Int_t i[]={4}
48 : /// AliTPCTransform trafo
49 : /// trafo.Transform(x,i,0,1)
50 : /// ~~~
51 : ///
52 : /// \author Magnus Mager, Marian Ivanov
53 :
54 : #include "AliTPCROC.h"
55 : #include "AliTPCCalPad.h"
56 : #include "AliTPCCalROC.h"
57 : #include "AliTPCcalibDB.h"
58 : #include "AliTPCParam.h"
59 : #include "TMath.h"
60 : #include "AliLog.h"
61 : #include "AliTPCExB.h"
62 : #include "AliTPCCorrection.h"
63 : #include "TGeoMatrix.h"
64 : #include "AliTPCRecoParam.h"
65 : #include "AliTPCCalibVdrift.h"
66 : #include "AliTPCTransform.h"
67 : #include "AliMagF.h"
68 : #include "TGeoGlobalMagField.h"
69 : #include "AliTracker.h"
70 : #include <AliCTPTimeParams.h>
71 : #include "AliTPCChebCorr.h"
72 : #include "AliTPCChebDist.h"
73 : #include "AliCDBManager.h"
74 : #include "AliCDBEntry.h"
75 : #include "TTreeStream.h"
76 : #include "AliLumiTools.h"
77 : #include "AliTPCclusterMI.h"
78 : #include <TGraph.h>
79 :
80 : /// \cond CLASSIMP
81 24 : ClassImp(AliTPCTransform)
82 : /// \endcond
83 :
84 :
85 24 : const Double_t AliTPCTransform::fgkSin20 = TMath::Sin(TMath::Pi()/9); // sin(20)
86 24 : const Double_t AliTPCTransform::fgkCos20 = TMath::Cos(TMath::Pi()/9); // cos(20)
87 24 : const Double_t AliTPCTransform::fgkMaxY2X = TMath::Tan(TMath::Pi()/18); // tg(10)
88 :
89 :
90 :
91 : AliTPCTransform::AliTPCTransform():
92 3 : AliTransform(),
93 3 : fCurrentRecoParam(0), //! current reconstruction parameters
94 3 : fCorrMapCacheRef(0),
95 3 : fCorrMapCache0(0),
96 3 : fCorrMapCache1(0),
97 3 : fCurrentMapScaling(1.0),
98 3 : fCorrMapLumiCOG(0.0),
99 3 : fLumiGraphRun(0),
100 3 : fLumiGraphMap(0),
101 3 : fCurrentRun(0), //! current run
102 3 : fCurrentTimeStamp(0), //! current time stamp
103 3 : fTimeDependentUpdated(kFALSE),
104 3 : fDebugStreamer(0)
105 15 : {
106 : //
107 : // Speed it up a bit!
108 : //
109 114 : for (Int_t i=0;i<18;++i) {
110 54 : Double_t alpha=TMath::DegToRad()*(10.+20.*(i%18));
111 54 : fSins[i]=TMath::Sin(alpha);
112 54 : fCoss[i]=TMath::Cos(alpha);
113 : }
114 3 : fPrimVtx[0]=0;
115 3 : fPrimVtx[1]=0;
116 3 : fPrimVtx[2]=0;
117 30 : for (int i=0;i<4;i++) {fLastCorr[i] = fLastCorrRef[i] = 0.0f;}
118 6 : }
119 :
120 : AliTPCTransform::AliTPCTransform(const AliTPCTransform& transform):
121 0 : AliTransform(transform),
122 0 : fCurrentRecoParam(transform.fCurrentRecoParam), //! current reconstruction parameters
123 0 : fCorrMapCacheRef(transform.fCorrMapCacheRef),
124 0 : fCorrMapCache0(transform.fCorrMapCache0),
125 0 : fCorrMapCache1(transform.fCorrMapCache1),
126 0 : fLumiGraphRun(transform.fLumiGraphRun ? new TGraph(*transform.fLumiGraphRun) : 0),
127 0 : fLumiGraphMap(transform.fLumiGraphMap ? new TGraph(*transform.fLumiGraphMap) : 0),
128 0 : fCurrentRun(transform.fCurrentRun), //! current run
129 0 : fCurrentTimeStamp(transform.fCurrentTimeStamp), //! current time stamp
130 0 : fTimeDependentUpdated(transform.fTimeDependentUpdated),
131 0 : fDebugStreamer(0)
132 0 : {
133 : /// Speed it up a bit!
134 :
135 0 : for (Int_t i=0;i<18;++i) {
136 0 : Double_t alpha=TMath::DegToRad()*(10.+20.*(i%18));
137 0 : fSins[i]=TMath::Sin(alpha);
138 0 : fCoss[i]=TMath::Cos(alpha);
139 : }
140 0 : fPrimVtx[0]=0;
141 0 : fPrimVtx[1]=0;
142 0 : fPrimVtx[2]=0;
143 0 : for (int i=0;i<4;i++) {fLastCorr[i] = fLastCorrRef[i] = 0.0f;}
144 0 : }
145 :
146 0 : AliTPCTransform::~AliTPCTransform() {
147 : /// Destructor
148 0 : delete fLumiGraphRun; // own copy should be attached
149 0 : delete fLumiGraphMap; // own copy should be attached
150 0 : delete fCorrMapCacheRef;
151 0 : delete fCorrMapCache0;
152 0 : delete fCorrMapCache1;
153 0 : }
154 :
155 : void AliTPCTransform::SetPrimVertex(Double_t *vtx){
156 : ///
157 :
158 0 : fPrimVtx[0]=vtx[0];
159 0 : fPrimVtx[1]=vtx[1];
160 0 : fPrimVtx[2]=vtx[2];
161 0 : }
162 :
163 :
164 : void AliTPCTransform::Transform(Double_t *x,Int_t *i,UInt_t /*time*/,
165 : Int_t /*coordinateType*/) {
166 : /// input: x[0] - pad row
167 : /// x[1] - pad
168 : /// x[2] - time in us
169 : /// i[0] - sector
170 : /// output: x[0] - x (all in the rotated global coordinate frame)
171 : /// x[1] - y
172 : /// x[2] - z
173 : ///
174 : /// primvtx - position of the primary vertex
175 : /// used for the TOF correction
176 : /// TOF of particle calculated assuming the speed-of-light and
177 : /// line approximation
178 :
179 259792 : if (!fCurrentRecoParam) return;
180 129896 : Int_t row=TMath::Nint(x[0]);
181 129896 : Int_t pad=TMath::Nint(x[1]);
182 129896 : Int_t sector=i[0];
183 129896 : AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
184 : //
185 129896 : AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
186 129896 : Double_t bzField = magF->SolenoidField(); //field in kGaus
187 :
188 129896 : AliTPCCalPad * time0TPC = calib->GetPadTime0();
189 129896 : AliTPCParam * param = calib->GetParameters();
190 :
191 129896 : if (!param){
192 0 : AliFatal("Parameters missing");
193 0 : return; // make coverity happy
194 : }
195 129896 : if (!time0TPC){
196 0 : AliFatal("Time unisochronity missing");
197 0 : return ; // make coverity happy
198 : }
199 :
200 : // Apply Time0 correction - Pad by pad fluctuation
201 259792 : if (!calib->HasAlignmentOCDB()) x[2]-=time0TPC->GetCalROC(sector)->GetValue(row,pad);
202 : //
203 : // Tranform from pad - time coordinate system to the rotated global (tracking) system
204 129896 : Local2RotatedGlobal(sector,x);
205 : Bool_t isInRotated = kTRUE;
206 : //
207 : //
208 129896 : if (fCurrentRecoParam->GetUseCorrectionMap()) ApplyCorrectionMap(sector, row, x);
209 : // Alignment
210 : //TODO: calib->GetParameters()->GetClusterMatrix(sector)->LocalToMaster(x,xx);
211 : //
212 129896 : if(fCurrentRecoParam->GetUseExBCorrection()) {
213 129896 : RotatedGlobal2Global(sector,x);
214 : isInRotated = kFALSE;
215 129896 : Double_t xx[3];
216 129896 : calib->GetExB()->Correct(x,xx); // old ExB correction
217 1039168 : for (int i=3;i--;) x[i] = xx[i];
218 129896 : }
219 :
220 : //
221 : // new composed correction - will replace soon ExB correction
222 : //
223 129896 : if(fCurrentRecoParam->GetUseComposedCorrection()) {
224 : //
225 0 : if (isInRotated) {
226 0 : RotatedGlobal2Global(sector,x);
227 : isInRotated = kFALSE;
228 0 : }
229 0 : AliTPCCorrection * correction = calib->GetTPCComposedCorrection(); // first user defined correction // if does not exist try to get it from calibDB array
230 0 : if (!correction) correction = calib->GetTPCComposedCorrection(bzField);
231 0 : AliTPCCorrection * correctionDelta = calib->GetTPCComposedCorrectionDelta();
232 0 : if (correction) {
233 0 : Float_t distPoint[3]={static_cast<Float_t>(x[0]),static_cast<Float_t>(x[1]),static_cast<Float_t>(x[2])};
234 0 : correction->CorrectPoint(distPoint, sector);
235 0 : for (int i=3;i--;) x[i]=distPoint[i];
236 0 : if (correctionDelta&&fCurrentRecoParam->GetUseAlignmentTime()){ // appply time dependent correction if available and enabled
237 0 : Float_t distPointDelta[3]={static_cast<Float_t>(x[0]),static_cast<Float_t>(x[1]),static_cast<Float_t>(x[2])};
238 0 : correctionDelta->CorrectPoint(distPointDelta, sector);
239 0 : for (int i=3;i--;) x[i]=distPointDelta[i];
240 0 : }
241 0 : }
242 0 : }
243 :
244 : //
245 : // Time of flight correction
246 : //
247 129896 : if (fCurrentRecoParam->GetUseTOFCorrection()){
248 129896 : const Int_t kNIS=param->GetNInnerSector(), kNOS=param->GetNOuterSector();
249 : Float_t sign=1;
250 129896 : if (sector < kNIS) {
251 92200 : sign = (sector < kNIS/2) ? 1 : -1;
252 92200 : } else {
253 37696 : sign = ((sector-kNIS) < kNOS/2) ? 1 : -1;
254 : }
255 : Float_t deltaDr =0;
256 : Float_t dist=0;
257 129896 : dist+=x[0]*x[0]; // (fPrimVtx[0]-x[0])*(fPrimVtx[0]-x[0]); // RS the vertex is anyway close to 0
258 129896 : dist+=x[1]*x[1]; // (fPrimVtx[1]-x[1])*(fPrimVtx[1]-x[1]);
259 129896 : dist+=(fPrimVtx[2]-x[2])*(fPrimVtx[2]-x[2]);
260 129896 : dist = TMath::Sqrt(dist);
261 : // drift length correction because of TOF
262 : // the drift velocity is in cm/s therefore multiplication by 0.01
263 129896 : deltaDr = (dist*(0.01*param->GetDriftV()))/TMath::C();
264 129896 : x[2]+=sign*deltaDr;
265 129896 : }
266 : //
267 : //
268 129896 : if (!isInRotated) {
269 129896 : Global2RotatedGlobal(sector,x);
270 : isInRotated = kTRUE;
271 129896 : }
272 : //
273 : // Apply non linear distortion correction
274 : //
275 129896 : int useFieldCorr = fCurrentRecoParam->GetUseFieldCorrection();
276 129896 : if (useFieldCorr&(0x2|0x4|0x8)) {
277 : AliTPCCalPad * distortionMapY=0,*distortionMapZ=0,*distortionMapR=0;
278 : //
279 : // wt - to get it form the OCDB
280 : // ignore T1 and T2
281 0 : Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
282 : Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
283 0 : if (sector%36<18) ezField*=-1;
284 : // Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
285 0 : Double_t wt = -10.0 * (bzField) * vdrift / ezField ; // RS: field is already in kGauss
286 0 : Double_t c0=1./(1.+wt*wt);
287 0 : Double_t c1=wt/c0;
288 :
289 : //can be switch on for each dimension separatelly
290 0 : if (useFieldCorr&0x2 && (distortionMapY=calib->GetDistortionMap(0))) {
291 0 : x[1]-= c0*distortionMapY->GetCalROC(sector)->GetValue(row,pad);
292 0 : x[0]-= c1*distortionMapY->GetCalROC(sector)->GetValue(row,pad);
293 0 : }
294 0 : if (useFieldCorr&0x4 && (distortionMapZ=calib->GetDistortionMap(1))) {
295 0 : x[2]-=distortionMapZ->GetCalROC(sector)->GetValue(row,pad);
296 0 : }
297 0 : if (useFieldCorr&0x8 && (distortionMapR=calib->GetDistortionMap(2))) {
298 0 : x[0]-= c0*distortionMapR->GetCalROC(sector)->GetValue(row,pad);
299 0 : x[1]-=-c1*distortionMapR->GetCalROC(sector)->GetValue(row,pad)*wt;
300 0 : }
301 : //
302 0 : }
303 : //
304 259792 : }
305 :
306 : void AliTPCTransform::Local2RotatedGlobal(Int_t sector, Double_t *x) const {
307 : /// Tranform coordinate from
308 : /// row, pad, time to x,y,z
309 : ///
310 : /// Drift Velocity
311 : /// Current implementation - common drift velocity - for full chamber
312 : /// TODO: use a map or parametrisation!
313 :
314 259792 : if (!fCurrentRecoParam) return;
315 : const Int_t kMax =60; // cache for 60 seconds
316 : static time_t lastStamp=-1; //cached values
317 : static Double_t lastCorr = 1;
318 : //
319 129896 : AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
320 129896 : AliTPCParam * param = calib->GetParameters();
321 129896 : AliTPCCalibVdrift *driftCalib = AliTPCcalibDB::Instance()->GetVdrift(fCurrentRun);
322 129896 : Double_t driftCorr = 1.;
323 129896 : if (driftCalib){
324 : //
325 : // caching drift correction - temp. fix
326 : // Extremally slow procedure
327 0 : if ( TMath::Abs(Int_t(lastStamp)-Int_t(fCurrentTimeStamp))<kMax){
328 0 : driftCorr = lastCorr;
329 0 : }else{
330 0 : driftCorr = 1.+(driftCalib->GetPTRelative(fCurrentTimeStamp,0)+ driftCalib->GetPTRelative(fCurrentTimeStamp,1))*0.5;
331 0 : lastCorr=driftCorr;
332 0 : lastStamp=fCurrentTimeStamp;
333 :
334 : }
335 : }
336 : //
337 : // simple caching non thread save
338 : static Double_t vdcorrectionTime=1;
339 : static Double_t vdcorrectionTimeGY=0;
340 : static Double_t time0corrTime=0;
341 : static Double_t deltaZcorrTime=0;
342 : static time_t lastStampT=-1;
343 : //
344 129896 : Bool_t isChange=(lastStampT!=(Int_t)fCurrentTimeStamp)||lastStamp<0;
345 129896 : if (lastStampT!=(Int_t)fCurrentTimeStamp){
346 2 : lastStampT=fCurrentTimeStamp;
347 2 : if(fCurrentRecoParam->GetUseDriftCorrectionTime()>0) {
348 4 : vdcorrectionTime = (1+AliTPCcalibDB::Instance()->
349 2 : GetVDriftCorrectionTime(fCurrentTimeStamp,
350 2 : fCurrentRun,
351 2 : sector%36>=18,
352 2 : fCurrentRecoParam->GetUseDriftCorrectionTime()));
353 4 : time0corrTime= AliTPCcalibDB::Instance()->
354 2 : GetTime0CorrectionTime(fCurrentTimeStamp,
355 2 : fCurrentRun,
356 : sector%36>=18,
357 2 : fCurrentRecoParam->GetUseDriftCorrectionTime());
358 : //
359 4 : deltaZcorrTime= AliTPCcalibDB::Instance()->
360 2 : GetVDriftCorrectionDeltaZ(fCurrentTimeStamp,
361 2 : fCurrentRun,
362 : sector%36>=18,
363 : 0);
364 :
365 2 : }
366 : //
367 2 : if(fCurrentRecoParam->GetUseDriftCorrectionGY()>0) {
368 :
369 4 : Double_t corrGy= AliTPCcalibDB::Instance()->
370 2 : GetVDriftCorrectionGy(fCurrentTimeStamp,
371 2 : AliTPCcalibDB::Instance()->GetRun(),
372 2 : sector%36>=18,
373 2 : fCurrentRecoParam->GetUseDriftCorrectionGY());
374 2 : vdcorrectionTimeGY = corrGy;
375 2 : }
376 : }
377 :
378 :
379 129896 : if (!param){
380 0 : AliFatal("Parameters missing");
381 0 : return; // make coverity happy
382 : }
383 129896 : Int_t row=TMath::Nint(x[0]);
384 : // Int_t pad=TMath::Nint(x[1]);
385 : //
386 129896 : const Int_t kNIS=param->GetNInnerSector(), kNOS=param->GetNOuterSector();
387 : Double_t sign = 1.;
388 129896 : Double_t zwidth = param->GetZWidth()*driftCorr;
389 129896 : Float_t xyzPad[3];
390 129896 : AliTPCROC::Instance()->GetPositionGlobal(sector, TMath::Nint(x[0]) ,TMath::Nint(x[1]), xyzPad);
391 259792 : if (AliTPCRecoParam:: GetUseTimeCalibration()) zwidth*=vdcorrectionTime*(1+xyzPad[1]*vdcorrectionTimeGY);
392 : Double_t padWidth = 0;
393 : Double_t padLength = 0;
394 : Double_t maxPad = 0;
395 : //
396 129896 : if (sector < kNIS) {
397 92200 : maxPad = param->GetNPadsLow(row);
398 92200 : sign = (sector < kNIS/2) ? 1 : -1;
399 92200 : padLength = param->GetPadPitchLength(sector,row);
400 92200 : padWidth = param->GetPadPitchWidth(sector);
401 92200 : } else {
402 37696 : maxPad = param->GetNPadsUp(row);
403 37696 : sign = ((sector-kNIS) < kNOS/2) ? 1 : -1;
404 37696 : padLength = param->GetPadPitchLength(sector,row);
405 37696 : padWidth = param->GetPadPitchWidth(sector);
406 : }
407 : //
408 : // X coordinate
409 129896 : x[0] = param->GetPadRowRadii(sector,row); // padrow X position - ideal
410 : //
411 : // Y coordinate
412 : //
413 129896 : x[1]=(x[1]-0.5*maxPad)*padWidth;
414 : // pads are mirrorred on C-side
415 129896 : if (sector%36>17){
416 62556 : x[1]*=-1;
417 62556 : }
418 :
419 : //
420 :
421 : //
422 : // Z coordinate
423 : //
424 129896 : Double_t delay=0;
425 129896 : if (AliTPCcalibDB::Instance()->IsTrgL0()){
426 : // by defualt we assume L1 trigger is used - make a correction in case of L0
427 0 : AliCTPTimeParams* ctp = AliTPCcalibDB::Instance()->GetCTPTimeParams();
428 0 : if (ctp){
429 : //for TPC standalone runs no ctp info
430 0 : delay = ctp->GetDelayL1L0()*0.000000025;
431 0 : delay/=param->GetTSample();
432 0 : x[2]-=delay;
433 0 : }
434 0 : }
435 259792 : if (isChange && fDebugStreamer!=NULL){
436 : //
437 : //
438 0 : Double_t zwidth=param->GetZWidth();
439 0 : Double_t binsL1=param->GetNTBinsL1();
440 0 : Double_t zsigma=param->GetZSigma();
441 0 : (*fDebugStreamer)<<"transformDump"<<
442 0 : "fCurrentTimeStamp="<<lastStampT<<
443 0 : "zwidth="<<zwidth<< // nominal z drift
444 0 : "delay="<<delay<< // trigger delay
445 0 : "binsL1="<<binsL1<< // L1 delay
446 0 : "zsigma="<<zsigma<< // z sigma
447 0 : "driftCorr="<<driftCorr<<
448 0 : "vdcorrectionTime="<<vdcorrectionTime<<
449 0 : "time0corrTime="<<time0corrTime<<
450 0 : "deltaZcorrTime="<<deltaZcorrTime<<
451 0 : "vdcorrectionTimeGY="<<vdcorrectionTimeGY<<
452 : "\n";
453 :
454 :
455 0 : }
456 129896 : x[2]-= param->GetNTBinsL1();
457 129896 : x[2]*= zwidth; // tranform time bin to the distance to the ROC
458 129896 : x[2]-= 3.*param->GetZSigma() + time0corrTime;
459 : // subtract the time offsets
460 129896 : x[2] = sign*( param->GetZLength(sector) - x[2]);
461 129896 : x[2]-=deltaZcorrTime; // subtrack time dependent z shift (calibrated together with the drift velocity and T0)
462 389688 : }
463 :
464 : void AliTPCTransform::RotatedGlobal2Global(Int_t sector,Double_t *x) const {
465 : /// transform possition rotated global to the global
466 :
467 259792 : Double_t cos,sin;
468 129896 : GetCosAndSin(sector,cos,sin);
469 129896 : Double_t tmp=x[0];
470 129896 : x[0]= cos*tmp-sin*x[1];
471 129896 : x[1]=+sin*tmp+cos*x[1];
472 129896 : }
473 :
474 : void AliTPCTransform::Global2RotatedGlobal(Int_t sector,Double_t *x) const {
475 : /// tranform possition Global2RotatedGlobal
476 :
477 259792 : Double_t cos,sin;
478 129896 : GetCosAndSin(sector,cos,sin);
479 129896 : Double_t tmp=x[0];
480 129896 : x[0]= cos*tmp+sin*x[1];
481 129896 : x[1]= -sin*tmp+cos*x[1];
482 129896 : }
483 :
484 : void AliTPCTransform::GetCosAndSin(Int_t sector,Double_t &cos,
485 : Double_t &sin) const {
486 519584 : cos=fCoss[sector%18];
487 259792 : sin=fSins[sector%18];
488 259792 : }
489 :
490 :
491 : void AliTPCTransform::ApplyTransformations(Double_t */*xyz*/, Int_t /*volID*/){
492 : /// Modify global position
493 : /// xyz - global xyz position
494 : /// volID - volID of detector (sector number)
495 :
496 0 : }
497 :
498 : void AliTPCTransform::SetCurrentTimeStamp(time_t timeStamp)
499 : {
500 : // set event time stamp and if needed, upload caches
501 42 : fCurrentTimeStamp = timeStamp;
502 21 : fTimeDependentUpdated = kFALSE;
503 21 : UpdateTimeDependentCache();
504 21 : }
505 :
506 : Bool_t AliTPCTransform::UpdateTimeDependentCache()
507 : {
508 : // update cache for time-dependent parameters
509 : //
510 : static time_t lastTimeStamp = -1;
511 42 : fTimeDependentUpdated = kFALSE;
512 : //
513 21 : Bool_t timeChanged = lastTimeStamp!=fCurrentTimeStamp;
514 21 : if (!fCurrentRecoParam) {
515 0 : AliWarning("RecoParam is not set, do nothing");
516 0 : return fTimeDependentUpdated;
517 : }
518 21 : while (fCurrentRecoParam->GetUseCorrectionMap()) {
519 0 : if (!fCorrMapCacheRef) fCorrMapCacheRef = LoadFieldDependendStaticCorrectionMap(kTRUE); // need to load the reference correction map
520 : //
521 0 : int mapTimeDepMethod = fCurrentRecoParam->GetCorrMapTimeDepMethod();
522 0 : Bool_t needToLoad = timeChanged;
523 0 : if (fCorrMapCache0 && timeChanged) { // do we need to update already loaded maps?
524 : // 1: easiest case: map is already cached, it is either time-static or there is
525 : // no other map to follow
526 0 : if (!fCorrMapCache0->GetTimeDependent()) needToLoad = kFALSE; // maps are not time dependent, no update needed
527 : //
528 : // 2: still easy: time dependent in interpolation mode but we have in memory all what we need
529 0 : if (fCorrMapCache1 && // interpolation method already triggered loading of 2nd map
530 0 : ( (fCurrentTimeStamp>=fCorrMapCache0->GetTimeStampCenter()&& // still covered by already loaded
531 0 : fCurrentTimeStamp<=fCorrMapCache1->GetTimeStampCenter()) // interpolation region
532 0 : ||
533 0 : (fCurrentTimeStamp<fCorrMapCache0->GetTimeStampEnd()&&fCorrMapCache0->IsFirst()) // no earlier object
534 0 : ||
535 0 : (fCurrentTimeStamp>fCorrMapCache1->GetTimeStampStart()&&fCorrMapCache1->IsLast()) // no later object
536 : ) ) {
537 : needToLoad = kFALSE; // no update needed
538 0 : }
539 : // 3: still easy: time depenedent but only 1 map needed (no interpolation)
540 0 : if ( (mapTimeDepMethod!=AliTPCRecoParam::kCorrMapInterpolation)
541 0 : &&
542 0 : ((fCurrentTimeStamp<=fCorrMapCache0->GetTimeStampEnd() &&
543 0 : (fCurrentTimeStamp>=fCorrMapCache0->GetTimeStampStart()||fCorrMapCache0->IsFirst()))
544 0 : ||
545 0 : (fCurrentTimeStamp>=fCorrMapCache0->GetTimeStampStart() &&
546 0 : (fCurrentTimeStamp<=fCorrMapCache0->GetTimeStampEnd()||fCorrMapCache0->IsLast()))
547 : )) {
548 : needToLoad = kFALSE; // still covered by existing map or no other map to load
549 0 : }
550 : }
551 : //
552 0 : if (needToLoad) { // need to upload correction maps, potentially time dependent
553 0 : TObjArray* mapsArr = LoadCorrectionMaps(kFALSE);
554 : // are these time-static maps?
555 0 : if (!((AliTPCChebCorr*)mapsArr->UncheckedAt(0))->GetTimeDependent()) {
556 0 : fCorrMapCache0 = LoadFieldDependendStaticCorrectionMap(kFALSE,mapsArr); // static maps are field-dependent
557 0 : }
558 : else {
559 0 : LoadCorrectionMapsForTimeBin(mapsArr); // load maps matching to time stamp
560 : }
561 : // clean unneeded object to save the memory
562 0 : mapsArr->Remove((TObject*)fCorrMapCache0);
563 0 : if (fCorrMapCache1) mapsArr->Remove((TObject*)fCorrMapCache1);
564 0 : mapsArr->SetOwner(kTRUE);
565 0 : delete mapsArr;
566 : //
567 0 : if (fCorrMapCache0 && !fCorrMapCache0->IsCorrection())
568 0 : AliFatalF("Uploaded map is not correction: %s",fCorrMapCache0->IsA()->GetName());
569 0 : if (fCorrMapCache1 && !fCorrMapCache1->IsCorrection())
570 0 : AliFatalF("Uploaded map is not correction: %s",fCorrMapCache1->IsA()->GetName());
571 :
572 : // check time stamps
573 0 : if (fCorrMapCache0 && fCorrMapCache0->GetTimeStampStart()>fCurrentTimeStamp) {
574 0 : AliWarningF("Event timestamp %ld < map0 beginning %ld",fCurrentTimeStamp,fCorrMapCache0->GetTimeStampStart());
575 0 : }
576 0 : if (fCorrMapCache1) {
577 0 : if (fCorrMapCache1->GetTimeStampEnd()<fCurrentTimeStamp) {
578 0 : AliWarningF("Event timestamp %ld > map1 end %ld",fCurrentTimeStamp,fCorrMapCache1->GetTimeStampEnd());
579 0 : }
580 : }
581 0 : else if (fCorrMapCache0->GetTimeStampEnd()<fCurrentTimeStamp) {
582 0 : AliWarningF("Event timestamp %ld > map0 end %ld",fCurrentTimeStamp,fCorrMapCache0->GetTimeStampEnd());
583 0 : }
584 0 : }
585 : // do we need to update luminosity scaling params?
586 0 : if (timeChanged) {
587 0 : switch(mapTimeDepMethod) {
588 : case AliTPCRecoParam::kCorrMapInterpolation :
589 : case AliTPCRecoParam::kCorrMapNoScaling : break;
590 : case AliTPCRecoParam::kCorrMapGlobalScalingLumi:
591 : // load lumi graph used for map (run may be different from current one if the map is default one)
592 0 : if (!fLumiGraphMap) {
593 : //if (!fLumiGraphMap || fLumiGraphMap->GetUniqueID()!=fCorrMapCache0->GetRun()) {
594 : //if (fLumiGraphMap) delete fLumiGraphMap;
595 0 : int runMap = fCorrMapCache0->GetRun();
596 0 : if (runMap<0) AliErrorF("CorrectionMap Lumi scaling requested but run is not set, current %d will be used",fCurrentRun);
597 0 : fLumiGraphMap = AliLumiTools::GetLumiGraph(fCurrentRecoParam->GetUseLumiType(),runMap);
598 0 : if (!fLumiGraphMap) AliFatalF("Failed to load CorrectionMapluminosity lumi graph of type %d",fCurrentRecoParam->GetUseLumiType());
599 0 : }
600 : // load lumi graph for this run
601 0 : if (!fLumiGraphRun) {
602 0 : fLumiGraphRun = AliLumiTools::GetLumiGraph(fCurrentRecoParam->GetUseLumiType(),fCurrentRun);
603 0 : if (!fLumiGraphRun) AliFatalF("Failed to load run lumi graph of type %d",fCurrentRecoParam->GetUseLumiType());
604 : }
605 0 : if (needToLoad) { // calculate average luminosity used for current corr. map
606 0 : fCorrMapLumiCOG = fCorrMapCache0->GetLuminosityCOG(fLumiGraphMap); // recalculate lumi for new map
607 0 : if (!fCorrMapLumiCOG) AliError("Correction map rescaling with luminosity cannot be done, will use constant map");
608 : }
609 : //
610 0 : fCurrentMapScaling = 1.0;
611 0 : if (fCorrMapLumiCOG>0) {
612 0 : double lumi = fLumiGraphRun->Eval(fCurrentTimeStamp);
613 0 : fCurrentMapScaling = lumi/fCorrMapLumiCOG;
614 0 : AliInfoF("Luminosity scaling factor %.3f will be used for time %ld (Lumi: current: %.2e COG:%.2e)",
615 : fCurrentMapScaling,fCurrentTimeStamp,lumi,fCorrMapLumiCOG);
616 0 : }
617 : //
618 : break;
619 0 : default: AliFatalF("Unknown corrections time-evolution mode %d",mapTimeDepMethod);
620 0 : };
621 : }
622 : //
623 : break;
624 : } // loading of correction maps
625 : //
626 : // other time dependent stuff if needed
627 : //
628 21 : lastTimeStamp = fCurrentTimeStamp;
629 21 : fTimeDependentUpdated = kTRUE;
630 21 : return fTimeDependentUpdated;
631 21 : }
632 :
633 : //______________________________________________________
634 : void AliTPCTransform::LoadCorrectionMapsForTimeBin(TObjArray* mapsArrProvided)
635 : {
636 : // loads time-independent correction map for given time bin
637 : //
638 : TObjArray* mapsArr = mapsArrProvided;
639 0 : if (!mapsArr) mapsArr = LoadCorrectionMaps(kFALSE);
640 0 : int entries = mapsArr->GetEntriesFast();
641 0 : delete fCorrMapCache0;
642 0 : delete fCorrMapCache1;
643 0 : fCorrMapCache0 = fCorrMapCache1 = 0;
644 : //1) choose map covering current time stamp.
645 : int selID=0;
646 0 : for (selID=0;selID<entries;selID++) { // maps are ordered in time
647 0 : fCorrMapCache0 = (AliTPCChebCorr*)mapsArr->UncheckedAt(selID);
648 0 : if (fCurrentTimeStamp > fCorrMapCache0->GetTimeStampEnd()) {
649 0 : if (selID==entries-1) break; // in case the maps end limit < timestamp close to EOR
650 : }
651 0 : else if (fCurrentTimeStamp < fCorrMapCache0->GetTimeStampStart()) {
652 0 : if (!selID) break; // in case the maps start limit > timestamp close to SOR
653 : }
654 : else break; // exact match
655 : }
656 : // this should not happen
657 0 : if (!fCorrMapCache0) AliFatalF("Did not find map matching to timestamp %ld",fCurrentTimeStamp);
658 : //
659 0 : if (!selID) fCorrMapCache0->SetFirst(); // flag if there are maps before/after
660 0 : if (selID==entries-1) fCorrMapCache0->SetLast();
661 : //
662 : // if we are in the interpolation mode, load 2nd closest map
663 0 : if (fCurrentRecoParam->GetCorrMapTimeDepMethod()==AliTPCRecoParam::kCorrMapInterpolation) {
664 0 : if (fCurrentTimeStamp<=fCorrMapCache0->GetTimeStampCenter()) { // load map before or closest, if any
665 0 : if (selID) fCorrMapCache1 = (AliTPCChebCorr*)mapsArr->UncheckedAt(--selID);
666 0 : else if (selID<entries-1) fCorrMapCache1 = (AliTPCChebCorr*)mapsArr->UncheckedAt(++selID);
667 : }
668 : else { // load map after or closest, if any
669 0 : if (selID<entries-1) fCorrMapCache1 = (AliTPCChebCorr*)mapsArr->UncheckedAt(++selID);
670 0 : else if (selID) fCorrMapCache1 = (AliTPCChebCorr*)mapsArr->UncheckedAt(--selID);
671 : }
672 : //
673 0 : if (fCorrMapCache1) {
674 0 : if (!selID) fCorrMapCache1->SetFirst(); // flag if there are maps before/after
675 0 : if (selID==entries-1) fCorrMapCache0->SetLast();
676 0 : if (fCorrMapCache1->GetTimeStampStart()<fCorrMapCache0->GetTimeStampStart()) {
677 0 : AliTPCChebCorr* tmp = fCorrMapCache1; // swap to order in time
678 0 : fCorrMapCache1 = fCorrMapCache0;
679 0 : fCorrMapCache0 = tmp;
680 0 : }
681 : }
682 : //
683 : }
684 : //
685 0 : AliInfoF("Loaded %d maps for time stamp %ld",fCorrMapCache1?2:1,fCurrentTimeStamp);
686 0 : fCorrMapCache0->Print();
687 0 : if (fCorrMapCache1) fCorrMapCache1->Print();
688 : //
689 0 : if (!mapsArrProvided) { // if loaded locally, clean unnecessary stuff
690 0 : mapsArr->Remove((TObject*)fCorrMapCache0);
691 0 : if (fCorrMapCache1) mapsArr->Remove((TObject*)fCorrMapCache1);
692 0 : mapsArr->SetOwner(kTRUE);
693 0 : delete mapsArr;
694 : }
695 0 : }
696 :
697 : //______________________________________________________
698 : AliTPCChebCorr* AliTPCTransform::LoadFieldDependendStaticCorrectionMap(Bool_t ref, TObjArray* mapsArrProvided)
699 : {
700 : // loads time-independent correction map for relevan field polarity. If ref is true, then the
701 : // reference map is loaded
702 : //
703 : const float kZeroField = 0.1;
704 : TObjArray* mapsArr = mapsArrProvided;
705 0 : if (!mapsArr) mapsArr = LoadCorrectionMaps(ref);
706 0 : int entries = mapsArr->GetEntriesFast();
707 0 : AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); // think on extracting field once only
708 0 : Double_t bzField = magF->SolenoidField(); //field in kGaus
709 : Char_t expectType = AliTPCChebCorr::kFieldZero;
710 0 : if (bzField<-kZeroField) expectType = AliTPCChebCorr::kFieldNeg;
711 0 : else if (bzField> kZeroField) expectType = AliTPCChebCorr::kFieldPos;
712 : //
713 : AliTPCChebCorr* cormap = 0;
714 0 : for (int i=0;i<entries;i++) {
715 0 : AliTPCChebCorr* map = (AliTPCChebCorr*)mapsArr->At(i); if (!map) continue;
716 0 : if (!map->IsCorrection()) continue;
717 0 : Char_t mtp = map->GetFieldType();
718 0 : if (mtp==expectType || mtp==AliTPCChebCorr::kFieldAny) cormap = map;
719 0 : if (mtp==expectType) break;
720 0 : }
721 0 : if (!cormap) AliFatalGeneralF("AliTPCTransform","Did not find %s correction map",ref ? "reference":"");
722 :
723 0 : AliInfoGeneralF("AliTPCTransform","Loaded %s correction map",ref ? "reference":"");
724 0 : cormap->Print();
725 0 : if (cormap->GetFieldType() == AliTPCChebCorr::kFieldAny) {
726 0 : AliWarningGeneralF("AliTPCTransform","ATTENTION: no map for field %+.1f was found, placeholder map is used",bzField);
727 0 : }
728 : //
729 0 : cormap->SetFirst(); // flag absence of maps before
730 0 : cormap->SetLast(); // and after
731 : //
732 0 : if (!mapsArrProvided) { // if loaded locally, clean unnecessary stuff
733 0 : mapsArr->Remove((TObject*)cormap);
734 0 : mapsArr->SetOwner(kTRUE);
735 0 : delete mapsArr;
736 : }
737 0 : return cormap;
738 0 : }
739 :
740 : //______________________________________________________
741 : TObjArray* AliTPCTransform::LoadCorrectionMaps(Bool_t refMap)
742 : {
743 : // TPC fast Chebyshev correction map, loaded on demand, not handler by calibDB
744 : const char* kNameRef = "TPC/Calib/CorrectionMapsRef";
745 : const char* kNameRun = "TPC/Calib/CorrectionMaps";
746 0 : const char* mapTypeName = refMap ? kNameRef : kNameRun;
747 : //
748 0 : AliCDBManager* man = AliCDBManager::Instance();
749 0 : AliCDBEntry* entry = man->Get(mapTypeName);
750 0 : TObjArray* correctionMaps = (TObjArray*)entry->GetObject();
751 0 : for (int i=correctionMaps->GetEntriesFast();i--;) {
752 0 : AliTPCChebCorr* map = (AliTPCChebCorr*)correctionMaps->At(i);
753 0 : map->Init();
754 : }
755 0 : entry->SetOwner(0);
756 0 : entry->SetObject(0);
757 0 : man->UnloadFromCache(mapTypeName);
758 0 : return correctionMaps;
759 : //
760 0 : }
761 :
762 : //______________________________________________________
763 : void AliTPCTransform::ApplyCorrectionMap(int roc, int row, double xyzSect[3])
764 : {
765 : // apply correction from the map to a point at given ROC and row (IROC/OROC convention)
766 0 : EvalCorrectionMap(roc, row, xyzSect, fLastCorrRef, kTRUE);
767 0 : EvalCorrectionMap(roc, row, xyzSect, fLastCorr, kFALSE);
768 0 : if (fLastCorr[3]<1e-6) { // run specific map had no parameterization for this region, override by default
769 0 : for (int i=3;i--;) fLastCorr[i] = fLastCorrRef[i];
770 0 : fLastCorr[3] = 0.f;
771 0 : }
772 : else {
773 0 : fLastCorr[3] = fLastCorr[3]>fLastCorrRef[3] ? TMath::Sqrt(fLastCorr[3]*fLastCorr[3] - fLastCorrRef[3]*fLastCorrRef[3]) : 0;
774 0 : if (fCurrentMapScaling!=1.0f) {
775 0 : for (int i=3;i--;) fLastCorr[i] = (fLastCorr[i]-fLastCorrRef[i])*fCurrentMapScaling + fLastCorrRef[i];
776 0 : fLastCorr[3] *= fCurrentMapScaling;
777 0 : }
778 : }
779 0 : for (int i=3;i--;) xyzSect[i] += fLastCorr[i];
780 : //
781 0 : }
782 :
783 : //______________________________________________________
784 : Float_t AliTPCTransform::GetCorrMapComponent(int roc, int row, const double xyz[3], int dimOut)
785 : {
786 : // calculate final correction component from the map
787 0 : float corr = EvalCorrectionMap(roc, row, xyz, dimOut, kFALSE); // run specific correction
788 : // do we need to rescale?
789 0 : if (fCurrentMapScaling!=1.0f) {
790 0 : float corrRef = EvalCorrectionMap(roc, row, xyz, dimOut, kTRUE); // ref correction
791 0 : if (dimOut<3) corr = (corr-corrRef)*fCurrentMapScaling + corrRef; // rescale with lumi
792 0 : }
793 0 : return corr;
794 : }
795 :
796 : //______________________________________________________
797 : void AliTPCTransform::EvalCorrectionMap(int roc, int row, const double xyz[3], float *res, Bool_t ref)
798 : {
799 : // get correction from the map for a point at given ROC and row (IROC/OROC convention)
800 0 : if (!fTimeDependentUpdated && !UpdateTimeDependentCache()) AliFatal("Failed to update time-dependent cache");
801 :
802 0 : AliTPCChebCorr* map = ref ? fCorrMapCacheRef : fCorrMapCache0;
803 0 : float y2x=xyz[1]/xyz[0], z2x = map->GetUseZ2R() ? xyz[2]/xyz[0] : xyz[2];
804 0 : map->Eval(roc,row,y2x,z2x,res);
805 :
806 0 : if (ref) return; // this was time-independent ref.map query request
807 : //
808 : // for time dependent correction need to evaluate 2 maps, assuming linear dependence
809 0 : if (fCorrMapCache1) {
810 0 : float delta1[4]={0};
811 0 : fCorrMapCache1->Eval(roc,row,y2x,z2x,delta1);
812 0 : UInt_t t0 = fCorrMapCache0->GetTimeStampCenter();
813 0 : UInt_t t1 = fCorrMapCache1->GetTimeStampCenter();
814 : // possible division by 0 is checked at upload of maps
815 0 : double dtScale = (fCurrentTimeStamp-t0)/double(t1-t0);
816 0 : for (int i=4;i--;) res[i] += (delta1[i]-res[i])*dtScale;
817 0 : }
818 : //
819 0 : }
820 :
821 : //______________________________________________________
822 : Float_t AliTPCTransform::EvalCorrectionMap(int roc, int row, const double xyz[3], int dimOut, Bool_t ref)
823 : {
824 : // get correction for dimOut-th dimension from the map for a point at given ROC and row (IROC/OROC convention)
825 0 : if (!fTimeDependentUpdated && !UpdateTimeDependentCache()) AliFatal("Failed to update time-dependent cache");
826 :
827 0 : AliTPCChebCorr* map = ref ? fCorrMapCacheRef : fCorrMapCache0;
828 0 : float y2x=xyz[1]/xyz[0], z2x = map->GetUseZ2R() ? xyz[2]/xyz[0] : xyz[2];
829 0 : float corr = map->Eval(roc,row,y2x,z2x,dimOut);
830 : //
831 0 : if (ref) return corr; // this was time-independent query request
832 : //
833 : // for time dependent correction need to evaluate 2 maps, assuming linear dependence
834 0 : if (fCorrMapCache1) {
835 0 : float corr1 = fCorrMapCache1->Eval(roc,row,y2x,z2x,dimOut);
836 0 : UInt_t t0 = fCorrMapCache0->GetTimeStampCenter();
837 0 : UInt_t t1 = fCorrMapCache1->GetTimeStampCenter();
838 : // possible division by 0 is checked at upload of maps
839 0 : double dtScale = (t1-fCurrentTimeStamp)/double(t1-t0);
840 0 : corr += (corr1-corr)*dtScale;
841 0 : }
842 : //
843 0 : return corr;
844 0 : }
845 :
846 : //______________________________________________________
847 : void AliTPCTransform::EvalDistortionMap(int roc, const double xyzSector[3], float *res)
848 : {
849 : // get distortions from the map for a point at given ROC
850 0 : if (!fTimeDependentUpdated && !UpdateTimeDependentCache()) AliFatal("Failed to update time-dependent cache");
851 0 : if (!fCorrMapCache0->IsDistortion()) AliFatalF("Uploaded map is not distortion: %s",fCorrMapCache0->IsA()->GetName());
852 0 : float y2x=xyzSector[1]/xyzSector[0], z2x = fCorrMapCache0->GetUseZ2R() ? xyzSector[2]/xyzSector[0] : xyzSector[2];
853 0 : ((AliTPCChebDist*)fCorrMapCache0)->Eval(roc,xyzSector[0],y2x,z2x,res);
854 : //
855 : // for time dependent correction need to evaluate 2 maps, assuming linear dependence
856 0 : if (fCorrMapCache1) {
857 0 : float delta1[4] = {0.0f};
858 0 : ((AliTPCChebDist*)fCorrMapCache1)->Eval(roc,xyzSector[0],y2x,z2x,res);
859 0 : UInt_t t0 = fCorrMapCache0->GetTimeStampCenter();
860 0 : UInt_t t1 = fCorrMapCache1->GetTimeStampCenter();
861 : // possible division by 0 is checked at upload of maps
862 0 : double dtScale = (t1-fCurrentTimeStamp)/double(t1-t0);
863 0 : for (int i=4;i--;) res[i] += (delta1[i]-res[i])*dtScale;
864 0 : }
865 : //
866 0 : }
867 :
868 : //______________________________________________________
869 : void AliTPCTransform::ApplyDistortionMap(int roc, double xyzLab[3])
870 : {
871 : // apply distortion from the map to a point provided in LAB coordinate
872 : // at given ROC and row (IROC/OROC convention)
873 : double xyzSect[3];
874 0 : float res[3];
875 0 : Global2RotatedGlobal(roc,xyzLab);
876 0 : EvalDistortionMap(roc, xyzLab, res); // now we are in sector coordinates
877 0 : for (int i=3;i--;) xyzLab[i] += res[i];
878 0 : RotatedGlobal2Global(roc,xyzLab);
879 : //
880 0 : }
881 :
882 : //_____________________________________________________________________________________
883 : void AliTPCTransform::ErrY2Z2Syst(const AliTPCclusterMI * cl, const double tgPhi, const double tgLam,
884 : double &serry2, double &serrz2)
885 : {
886 : // static function to calculate systematic errorY^2 on the cluster with current reco param
887 : // Must be static to be used also by external routines
888 : // tgPhi is the ab.tg(inclination wrt padrow), tgLam is tg of dip angle
889 : const float kEpsZBoundary = 1.e-6; // to disentangle A,C and A+C-common regions
890 : const float kMaxExpArg = 9.; // limit r-dumped error to this exp. argument
891 209276 : const float kMaxExpArgZ = TMath::Sqrt(2.*kMaxExpArg);
892 104638 : serry2 = serrz2 = 0;
893 : double sysErrY = 0, sysErrZ = 0;
894 : //
895 : // error on particular TPC subvolume (charging up)
896 104638 : const Double_t *errInner = fCurrentRecoParam->GetSystematicErrorClusterInner();
897 104638 : if (errInner[0]>0) {
898 : //
899 104638 : double dr = TMath::Abs(cl->GetX()-85.);
900 104638 : float argExp = dr/errInner[1]; // is the Z-range limited ?
901 104638 : if (argExp<kMaxExpArg) {
902 38342 : const TVectorF* zranges = fCurrentRecoParam->GetSystErrClInnerRegZ(); // is the Z-range of dumped error defined?
903 : int nz;
904 38342 : if (zranges && (nz=zranges->GetNoElements())) {
905 0 : Bool_t sideC = (cl->GetDetector()/18)&0x1; // is this C-side cluster
906 0 : const TVectorF* zrangesSigI = fCurrentRecoParam->GetSystErrClInnerRegZSigInv();
907 0 : float zcl = cl->GetZ(); // chose the closest range within 3 sigma
908 0 : const float* zr = zranges->GetMatrixArray();
909 0 : const float* zsi= zrangesSigI->GetMatrixArray();
910 : float dzarg = 999.;
911 0 : for (int i=nz;i--;) { // find closest region
912 0 : if (zr[i]<-kEpsZBoundary && !sideC) continue; // don't apply to A-side clusters if the Z-boundary is for C-region
913 0 : else if (zr[i]>kEpsZBoundary && sideC) continue; // don't apply to C-side clusters if the Z-boundary is for A-region
914 0 : float dz = TMath::Abs( (zr[i]-zcl)*zsi[i] );
915 0 : if (dz<dzarg) dzarg = dz;
916 : }
917 0 : if (dzarg<kMaxExpArgZ) { // is it small enough to call exp?
918 0 : argExp += 0.5*dzarg*dzarg;
919 0 : if (argExp<kMaxExpArg) sysErrY = sysErrZ = errInner[0]*TMath::Exp(-argExp);
920 : }
921 0 : }
922 38342 : else sysErrY = sysErrZ = errInner[0]*TMath::Exp(-argExp); // no condition on Z
923 38342 : }
924 :
925 104638 : }
926 104638 : serry2 += sysErrY*sysErrY;
927 104638 : serrz2 += sysErrZ*sysErrZ;
928 : //
929 418552 : sysErrY = fCurrentRecoParam->GetSystematicErrorClusterCustom() ?
930 209276 : fCurrentRecoParam->GetSystematicErrorClusterCustom()[0] : fCurrentRecoParam->GetSystematicErrorCluster()[0];
931 418552 : sysErrZ = fCurrentRecoParam->GetSystematicErrorClusterCustom() ?
932 209276 : fCurrentRecoParam->GetSystematicErrorClusterCustom()[1] : fCurrentRecoParam->GetSystematicErrorCluster()[1];
933 104638 : if (sysErrY) serry2 += sysErrY*sysErrY; // common syst error
934 104638 : if (sysErrZ) serrz2 += sysErrZ*sysErrZ; // common syst error
935 : //
936 104638 : double useDistY = fCurrentRecoParam->GetUseDistortionFractionAsErrorY();
937 104638 : if (useDistY>0) {
938 104638 : float dstY = cl->GetDistortionY()*useDistY;
939 104638 : float dstX = cl->GetDistortionX()*useDistY*tgPhi;
940 104638 : serry2 += dstY*dstY + dstX*dstX;
941 104638 : }
942 104638 : double useDispY = fCurrentRecoParam->GetUseDistDispFractionAsErrorY();
943 104638 : if (useDispY>0) {
944 104638 : useDispY *= cl->GetDistortionDispersion();
945 104638 : serry2 += useDispY*useDispY;
946 104638 : }
947 : //
948 104638 : double useDistZ = fCurrentRecoParam->GetUseDistortionFractionAsErrorZ();
949 104638 : if (useDistZ>0) {
950 104638 : float dstZ = cl->GetDistortionZ()*useDistZ;
951 104638 : float dstX = cl->GetDistortionX()*useDistZ*tgLam;
952 104638 : serrz2 += dstZ*dstZ + dstX*dstX;
953 104638 : }
954 104638 : double useDispZ = fCurrentRecoParam->GetUseDistDispFractionAsErrorZ();
955 104638 : if (useDispZ>0) {
956 104638 : useDispZ *= cl->GetDistortionDispersion();
957 104638 : serrz2 += useDispZ*useDispZ;
958 104638 : }
959 : //
960 104638 : }
961 :
962 : //_____________________________________________________________________________________
963 : Double_t AliTPCTransform::ErrY2Syst(const AliTPCclusterMI * cl, const double tgAng)
964 : {
965 : // static function to calculate systematic errorY^2 on the cluster with current reco param
966 : // Must be static to be used also by external routines
967 : // tgAng is the ab.tg(inclination wrt padrow)
968 : const float kEpsZBoundary = 1.e-6; // to disentangle A,C and A+C-common regions
969 : const float kMaxExpArg = 9.; // limit r-dumped error to this exp. argument
970 0 : const float kMaxExpArgZ = TMath::Sqrt(2.*kMaxExpArg);
971 : double sysErr2 = 0., sysErr = 0.;
972 : //
973 : // error on particular TPC subvolume (charging up)
974 0 : const Double_t *errInner = fCurrentRecoParam->GetSystematicErrorClusterInner();
975 0 : if (errInner[0]>0) {
976 : //
977 0 : double dr = TMath::Abs(cl->GetX()-85.);
978 0 : float argExp = dr/errInner[1]; // is the Z-range limited ?
979 0 : if (argExp<kMaxExpArg) {
980 0 : const TVectorF* zranges = fCurrentRecoParam->GetSystErrClInnerRegZ(); // is the Z-range of dumped error defined?
981 : int nz;
982 0 : if (zranges && (nz=zranges->GetNoElements())) {
983 0 : Bool_t sideC = (cl->GetDetector()/18)&0x1; // is this C-side cluster
984 0 : const TVectorF* zrangesSigI = fCurrentRecoParam->GetSystErrClInnerRegZSigInv();
985 0 : float zcl = cl->GetZ(); // chose the closest range within 3 sigma
986 0 : const float* zr = zranges->GetMatrixArray();
987 0 : const float* zsi= zrangesSigI->GetMatrixArray();
988 : float dzarg = 999.;
989 0 : for (int i=nz;i--;) { // find closest region
990 0 : if (zr[i]<-kEpsZBoundary && !sideC) continue; // don't apply to A-side clusters if the Z-boundary is for C-region
991 0 : else if (zr[i]>kEpsZBoundary && sideC) continue; // don't apply to C-side clusters if the Z-boundary is for A-region
992 0 : float dz = TMath::Abs( (zr[i]-zcl)*zsi[i] );
993 0 : if (dz<dzarg) dzarg = dz;
994 : }
995 0 : if (dzarg<kMaxExpArgZ) { // is it small enough to call exp?
996 0 : argExp += 0.5*dzarg*dzarg;
997 0 : if (argExp<kMaxExpArg) sysErr = errInner[0]*TMath::Exp(-argExp);
998 : }
999 0 : }
1000 0 : else sysErr = errInner[0]*TMath::Exp(-argExp); // no condition on Z
1001 0 : }
1002 :
1003 0 : }
1004 0 : sysErr2 += sysErr*sysErr;
1005 : //
1006 0 : sysErr = fCurrentRecoParam->GetSystematicErrorClusterCustom() ?
1007 0 : fCurrentRecoParam->GetSystematicErrorClusterCustom()[0] : fCurrentRecoParam->GetSystematicErrorCluster()[0];
1008 0 : if (sysErr) sysErr2 += sysErr*sysErr; // common syst error
1009 : //
1010 0 : double useDist = fCurrentRecoParam->GetUseDistortionFractionAsErrorY();
1011 0 : if (useDist>0) {
1012 0 : float dstY = cl->GetDistortionY()*useDist;
1013 0 : float dstX = cl->GetDistortionX()*useDist*tgAng;
1014 0 : sysErr2 += dstY*dstY + dstX*dstX;
1015 0 : }
1016 0 : double useDisp = fCurrentRecoParam->GetUseDistDispFractionAsErrorY();
1017 0 : if (useDisp>0) {
1018 0 : useDisp *= cl->GetDistortionDispersion();
1019 0 : sysErr2 += useDisp*useDisp;
1020 0 : }
1021 0 : return sysErr2;
1022 : }
1023 :
1024 : //_____________________________________________________________________________________
1025 : Double_t AliTPCTransform::ErrZ2Syst(const AliTPCclusterMI * cl, const double tgAng)
1026 : {
1027 : // static function to calculate systematic errorY^2 on the cluster with current reco param
1028 : // Must be static to be used also by external routines
1029 : // tgAng is the ab.tg(lambda)
1030 : const float kEpsZBoundary = 1.e-6; // to disentangle A,C and A+C-common regions
1031 : const float kMaxExpArg = 9.; // limit r-dumped error to this exp. argument
1032 0 : const float kMaxExpArgZ = TMath::Sqrt(2.*kMaxExpArg);
1033 : double sysErr2 = 0., sysErr = 0.;
1034 : //
1035 : // error on particular TPC subvolume (charging up)
1036 0 : const Double_t *errInner = fCurrentRecoParam->GetSystematicErrorClusterInner();
1037 0 : if (errInner[0]>0) {
1038 : //
1039 0 : double dr = TMath::Abs(cl->GetX()-85.);
1040 0 : float argExp = dr/errInner[1]; // is the Z-range limited ?
1041 0 : if (argExp<kMaxExpArg) {
1042 0 : const TVectorF* zranges = fCurrentRecoParam->GetSystErrClInnerRegZ(); // is the Z-range of dumped error defined?
1043 : int nz;
1044 0 : if (zranges && (nz=zranges->GetNoElements())) {
1045 0 : Bool_t sideC = (cl->GetDetector()/18)&0x1; // is this C-side cluster
1046 0 : const TVectorF* zrangesSigI = fCurrentRecoParam->GetSystErrClInnerRegZSigInv();
1047 0 : float zcl = cl->GetZ(); // chose the closest range within 3 sigma
1048 0 : const float* zr = zranges->GetMatrixArray();
1049 0 : const float* zsi= zrangesSigI->GetMatrixArray();
1050 : float dzarg = 999.;
1051 0 : for (int i=nz;i--;) { // find closest region
1052 0 : if (zr[i]<-kEpsZBoundary && !sideC) continue; // don't apply to A-side clusters if the Z-boundary is for C-region
1053 0 : else if (zr[i]>kEpsZBoundary && sideC) continue; // don't apply to C-side clusters if the Z-boundary is for A-region
1054 0 : float dz = TMath::Abs( (zr[i]-zcl)*zsi[i] );
1055 0 : if (dz<dzarg) dzarg = dz;
1056 : }
1057 0 : if (dzarg<kMaxExpArgZ) { // is it small enough to call exp?
1058 0 : argExp += 0.5*dzarg*dzarg;
1059 0 : if (argExp<kMaxExpArg) sysErr = errInner[0]*TMath::Exp(-argExp);
1060 : }
1061 0 : }
1062 0 : else sysErr = errInner[0]*TMath::Exp(-argExp); // no condition on Z
1063 0 : }
1064 :
1065 0 : }
1066 0 : sysErr2 += sysErr*sysErr;
1067 : //
1068 0 : sysErr = fCurrentRecoParam->GetSystematicErrorClusterCustom() ?
1069 0 : fCurrentRecoParam->GetSystematicErrorClusterCustom()[1] : fCurrentRecoParam->GetSystematicErrorCluster()[1];
1070 0 : if (sysErr) sysErr2 += sysErr*sysErr; // common syst error
1071 : //
1072 0 : double useDist = fCurrentRecoParam->GetUseDistortionFractionAsErrorZ();
1073 0 : if (useDist>0) {
1074 0 : float dstZ = cl->GetDistortionZ()*useDist;
1075 0 : float dstX = cl->GetDistortionX()*useDist*tgAng;
1076 0 : sysErr2 += dstZ*dstZ + dstX*dstX;
1077 0 : }
1078 0 : double useDisp = fCurrentRecoParam->GetUseDistDispFractionAsErrorZ();
1079 0 : if (useDisp>0) {
1080 0 : useDisp *= cl->GetDistortionDispersion();
1081 0 : sysErr2 += useDisp*useDisp;
1082 0 : }
1083 0 : return sysErr2;
1084 : }
1085 :
|
