Line data Source code
1 : #ifndef ALITPCCALIBRES_H
2 : #define ALITPCCALIBRES_H
3 : #include <TSystem.h>
4 : #include <TNamed.h>
5 : #include <TTree.h>
6 : #include <TBranch.h>
7 : #include <TFile.h>
8 : #include <TVectorF.h>
9 : #include <TVectorD.h>
10 : #include <TString.h>
11 : #include <TMath.h>
12 : #include <TGeoGlobalMagField.h>
13 : #include <TGrid.h>
14 : #include <TNDArray.h>
15 : #include <THn.h>
16 : #include <TH1F.h>
17 : #include <TF1.h>
18 : #include <TEnv.h>
19 : #include <TStopwatch.h>
20 : #include <TGraphErrors.h>
21 : #include <TGeoMatrix.h>
22 : #include <TH2F.h>
23 : #include <TProfile2D.h>
24 : #include "AliLog.h"
25 : #include "AliExternalTrackParam.h"
26 : #include "AliTPCcalibAlignInterpolation.h"
27 : #include "AliGeomManager.h"
28 : #include "AliCDBManager.h"
29 : #include "AliGRPManager.h"
30 : #include "AliMagF.h"
31 : #include "AliSysInfo.h"
32 : #include "TStatToolkit.h"
33 : #include "AliSymMatrix.h"
34 : #include "AliTPCChebCorr.h"
35 :
36 :
37 0 : class AliTPCDcalibRes: public TNamed
38 : {
39 : public:
40 : enum {kEpanechnikovKernel, kGaussianKernel}; // defined kernels
41 : enum {kNSect=18,kNSect2=2*kNSect,kNROC=4*kNSect,kNPadRows=159, kNRowIROC=63, kNRowOROC1=64, kNRowOROC2=32};
42 : enum {kAlignmentBugFixedBit = AliTPCcalibAlignInterpolation::kAlignmentBugFixedBit};
43 : enum {kVDriftCalibMode, kDistExtractMode, kDistClosureTestMode,kNCollectModes};
44 : enum {kDistDone=BIT(0),kDispDone=BIT(1),kSmoothDone=BIT(2),kMasked=BIT(7)};
45 : enum {kVDWithTOFBC=BIT(14),kVDWithoutTOFBC=BIT(15)};
46 : enum {kUseTRDonly,kUseTOFonly,kUseITSonly,kUseTRDorTOF,kNExtDetComb}; // which points to use
47 : enum {kSmtLinDim=4, kMaxSmtDim=7}; // max size of matrix for smoothing, for pol1 and pol2 options
48 : enum {kCtrITS,kCtrTRD,kCtrTOF,kCtrBC0,kCtrNtr,kCtrNbr}; // control branches for test stat
49 : enum {kLargeTimeDiff=1000*3600,kMaxOnStack=10000}; // misc consts
50 :
51 : // the voxels are defined in following space
52 : enum {kVoxZ, // Z/X sector coordinates
53 : kVoxF, // y/x in sector coordinates
54 : kVoxX, // sector X coordinate
55 : kVoxV, // variable within the voxel (delta, stat, etc): last dimension of all THn histos
56 : kVoxHDim, kVoxDim=kVoxHDim-1};
57 :
58 : enum {kResX,kResY,kResZ,kResD,kResDim,kResDimG=kResDim-1}; // output dimensions
59 : //
60 : struct dtv_t { // struct for vdrift calibration loca residual
61 : Char_t side; // +/- for A/C
62 : Double32_t dz; //[-40.,40.,15] // +-kMaxResidZVD
63 : Double32_t drift; //[-260.,260.,14] drift distance
64 : Double32_t ylab; //[-260.,260.,14] y coordinate
65 : Double32_t r; //[80.,250.,7] radius
66 : Int_t t; // time stamp
67 : //
68 0 : dtv_t() {memset(this,0,sizeof(dtv_t));}
69 : };
70 : //
71 : struct dts_t { // struct for basic local residual
72 : Double32_t dy; //[-20.,20.,15] // [-kMaxResid,kMaxResid,14]
73 : Double32_t dz; //[-20.,20.,15] // [-kMaxResid,kMaxResid,14]
74 : Double32_t tgSlp; //[-2,2,14] //[kMaxTgSlp,kMaxTgSlp,14]
75 : UChar_t bvox[kVoxDim]; // voxel bin info: VoxF,kVoxX,kVoxZ
76 : //
77 0 : dts_t() {memset(this,0,sizeof(dts_t));}
78 : };
79 :
80 : // structure for closure test residuals
81 : struct dtc_t
82 : {
83 : Int_t t; // time stamp
84 : Double32_t dyR; //[-15.,15.,14]
85 : Double32_t dzR; //[-15.,15.,14]
86 : Double32_t dyC; //[-15.,15.,14]
87 : Double32_t dzC; //[-15.,15.,14]
88 : Double32_t q2pt;//[-3,3,14]
89 : Double32_t tgLam;//[-2.,2.,14]
90 : Double32_t tgSlp;//[-2,2,14]
91 : Float_t x; //[80,160,14]
92 : Float_t y; //[-50,50,14]
93 : Float_t z; //[-250,250,14]
94 : UChar_t bvox[kVoxDim]; // voxel bin info: kVoxQ,kVoxF,kVoxX,kVoxZ
95 : //
96 0 : dtc_t() {memset(this,0,sizeof(dtc_t));}
97 : };
98 :
99 : struct bres_t {
100 : Float_t D[kResDim]; // values of extracted distortions
101 : Float_t E[kResDim]; // their errors
102 : Float_t DS[kResDim]; // smoothed residual
103 : Float_t DC[kResDim]; // Cheb parameterized residual
104 : Float_t EXYCorr; // correlation between extracted X and Y
105 : Float_t dYSigMAD; // MAD estimator of dY sigma (dispersion after slope removal)
106 : Float_t dZSigLTM; // Z sigma from unbinned LTM estimator
107 : Float_t stat[kVoxHDim]; // statistics: averages of each voxel dimension + entries
108 : UChar_t bvox[kVoxDim]; // voxel identifier, here the bvox[0] shows number of Q bins used for Y
109 : UChar_t bsec; // sector ID (0-35)
110 : UChar_t flags; // status flag
111 : //
112 0 : bres_t() {memset(this,0,sizeof(bres_t));}
113 : };
114 :
115 : struct delta_t { // structure to organized the input from delta trees
116 : TVectorF *vecDYTRD,*vecDZTRD,*vecDYITS,*vecDZITS,*vecDYTOF,*vecDZTOF,*vecZ,*vecR,*vecSec,*vecPhi;
117 : AliExternalTrackParam* param;
118 : Double_t tofBC;
119 : Int_t nPrimTracks;
120 : Int_t timeStamp;
121 : UShort_t npValid;
122 : Char_t trdOK,tofOK,itsOK;
123 : //
124 0 : delta_t() {memset(this,0,sizeof(delta_t));}
125 : };
126 :
127 : public:
128 :
129 : AliTPCDcalibRes(int run=0,Long64_t tmin=0,Long64_t tmax=9999999999,const char* resList=0);
130 : virtual ~AliTPCDcalibRes();
131 :
132 : void CalibrateVDrift();
133 : void FitDrift();
134 0 : TProfile2D* GetHistoVDTimeInt() const {return fHVDTimeInt;}
135 0 : TH2F* GetHistoVDTimeCorr() const {return fHVDTimeCorr;}
136 0 : Int_t GetDeltaTVDrift() const {return fDeltaTVD;}
137 0 : Int_t GetSigmaTVDrift() const {return fSigmaTVD;}
138 0 : void SetDeltaTVDrift(int v=120) {fDeltaTVD = v;}
139 0 : void SetSigmaTVDrift(int v=600) {fSigmaTVD = v;}
140 : //
141 : void ProcessFromDeltaTrees();
142 : void ProcessFromLocalBinnedTrees();
143 : void ReProcessFromResVoxTree(const char* resTreeFile, Bool_t backup=kTRUE);
144 : void Save(const char* name=0);
145 : void SaveFitDrift();
146 :
147 : TTree* InitDeltaFile(const char* name, Bool_t connect=kTRUE, const char* treeName="delta");
148 : Bool_t EstimateChunkStatistics();
149 : Bool_t CollectDataStatistics();
150 : Bool_t EnoughStatForVDrift(int tstamp=-1,float maxHolesFrac=0.05);
151 : Int_t ParseInputList();
152 : void CloseTreeFile(TTree* dtree);
153 : void Init();
154 : void CollectData(const int mode = kDistExtractMode);
155 : void FillLocalResidualsTrees();
156 : void FillCorrectedResiduals();
157 : void FillDriftResidualsTrees();
158 : void ClosureTest();
159 : void CreateLocalResidualsTrees(int mode);
160 : void CloseLocalResidualsTrees(int mode);
161 : void ProcessResiduals();
162 : void ReProcessResiduals();
163 : void ProcessDispersions();
164 : void ProcessSectorResiduals(int is);
165 : void ReProcessSectorResiduals(int is);
166 : void ProcessSectorDispersions(int is);
167 : void ProcessVoxelResiduals(int np, float* tg, float *dy, float *dz, bres_t& voxRes);
168 : void ProcessVoxelDispersions(int np, const float* tg, float *dy, bres_t& voxRes);
169 : Int_t ValidateVoxels(int isect);
170 : //
171 : void InitFieldGeom(Bool_t field=kTRUE,Bool_t geom=kTRUE);
172 : THnF* CreateVoxelStatHisto(int sect);
173 : void LoadVDrift();
174 : void MakeVDriftOCDB(TString targetOCDBstorage);
175 : Float_t GetDriftCorrection(float z, float x, float phi, int rocID);
176 : Float_t tgpXY(float x, float y, float q2p, float bz);
177 :
178 0 : TVectorD* GetVDriftParam() const {return (TVectorD*)fVDriftParam;}
179 0 : TGraphErrors* GetVDriftGraph() const {return (TGraphErrors*)fVDriftGraph;}
180 :
181 : void WriteStatHistos();
182 : void LoadStatHistos();
183 : void WriteResTree();
184 : Bool_t LoadDataFromResTree(const char* resTreeFile);
185 :
186 : void FixAlignmentBug(int sect, float q2pt, float bz, float& alp, float& x, float &z, float &deltaY, float &deltaZ);
187 :
188 : Bool_t ValidateTrack();
189 : Bool_t CompareToHelix(float *resHelixY, float *resHelixZ);
190 :
191 : int CheckResiduals(Bool_t* mask,float &rmsLongMA);
192 :
193 0 : const char* GetVoxResFileName() const {return Form("%sTree.root",kResOut);}
194 :
195 : //------------------------------------ misc. stat. methods
196 : static Float_t FitPoly1Robust(int np, float* x, float* y, float* res, float* err, float ltmCut);
197 : static void FitCircle(int np, const float* x, const float* y,
198 : double &xc, double &yc, double &r, float* dy=0);
199 : static void DiffToMA(int np, const float *y, const int winLR, float* diffMA);
200 : static int DiffToLocLine(int np, const float* x, const float *y, const int nVoisin, float *diffY);
201 : static int DiffToMedLine(int np, const float* x, const float *y, const int nVoisin, float *diffY);
202 : static float RoFunc(int np, const float* x, const float* y, float b, float &aa);
203 : static Float_t SelKthMin(int k, int np, float* arr);
204 : static void medFit(int np, const float* x, const float* y, float &a, float &b, float *err=0, float delI=0.f);
205 : static Int_t* LTMUnbinnedSig(int np, const float *arr, TVectorF ¶ms , Float_t sigTgt, Float_t minFrac=0.7, Bool_t sorted=kFALSE);
206 : static Float_t MAD2Sigma(int np, float* y);
207 : static Bool_t FitPoly2(const float* x,const float* y, const float* w, int np, float *res, float *err=0);
208 : static Bool_t FitPoly1(const float* x,const float* y, const float* w, int np, float *res, float *err=0);
209 : static Bool_t GetTruncNormMuSig(double a, double b, double &mean, double &sig);
210 : static void TruncNormMod(double a, double b, double mu0, double sig0, double &muCf, double &sigCf);
211 : static Double_t GetLogL(TH1F* histo, int bin0, int bin1, double &mu, double &sig, double &logL0);
212 :
213 : //------------------------------------
214 :
215 :
216 : Int_t Smooth0(int isect);
217 : Bool_t GetSmoothEstimate(int isect, float x, float p, float z, int which, float *res, float *deriv=0);
218 : void SetKernelType(int tp=kEpanechnikovKernel, float bwX=2.1, float bwP=2.1, float bwZ=1.7,
219 : float scX=1.f,float scP=1.f,float scZ=1.f);
220 : Bool_t GetSmooth1D(double xQuery,double valerr[2],int np,const double* x,const double* y,const double* err,
221 : double wKernel,int kType=kGaussianKernel,Bool_t usePoly2=kFALSE,
222 : Bool_t xIncreasing=kTRUE, Float_t maxD2Range=3.0) const;
223 0 : Bool_t GetSmoothPol2(int i) const {return fSmoothPol2[i];}
224 0 : void SetSmoothPol2(int i,Bool_t v=kTRUE) {fSmoothPol2[i] = v;}
225 : //
226 : void CreateCorrectionObject();
227 : void InitBinning();
228 : Int_t GetXBin(float x);
229 : Int_t GetRowID(float x);
230 :
231 : Bool_t FindVoxelBin(int sectID, float x, float y, float z, UChar_t bin[kVoxHDim],float voxVars[kVoxHDim]);
232 0 : TH1F* GetEventsRateHisto() const {return fEvRateH;}
233 0 : TH1F* GetTracksRateHisto() const {return fTracksRate;}
234 0 : TH1F* GetTOFBCHisto() const {return fTOFBCTestH;}
235 0 : void SetTracksRateHisto(TH1F* h) {fTracksRate = h;} // needed for hacks
236 : TH1F* ExtractTrackRate() const;
237 0 : Float_t GetValidFracXBin(int sect, int bin) const {return fValidFracXBin[sect][bin];}
238 0 : Int_t GetNSmoothingFailed(int sect) const {return fNSmoothingFailedBins[sect];}
239 0 : Int_t GetNTracksUsed() const {return fNTrSelTot;}
240 0 : Int_t GetNTracksWithOutliers() const {return fNTrSelTotWO;}
241 0 : Int_t GetMinTrackToUse() const {return fMinTracksToUse;}
242 0 : Int_t GetMinTracksPerVDBin() const {return fMinTracksPerVDBin;}
243 0 : Int_t GetNTestTracks() const {return fNTestTracks;}
244 0 : Float_t GetTestStat(int row, int col) const {return fTestStat[row][col];}
245 0 : TGraph* GetLumiGraph() const {return (TGraph*) fLumiCTPGraph;}
246 0 : Float_t GetLumiCOG() const {return fLumiCOG;}
247 : //
248 : Int_t GetXBinExact(float x);
249 : Float_t GetY2X(int ix, int iy);
250 : Float_t GetY2XLow(int ix, int iy);
251 : Float_t GetDY2X(int ix);
252 : Float_t GetDY2XI(int ix);
253 : Float_t GetX(int i);
254 : Float_t GetXLow(int i);
255 : Float_t GetDX(int i);
256 : Float_t GetDXI(int i);
257 : Int_t GetY2XBinExact(float y2x, int ix);
258 : Int_t GetY2XBin(float y2x, int ix);
259 : Int_t GetZ2XBinExact(float z2x);
260 : Int_t GetZ2XBin(float z2x);
261 : Float_t GetZ2XLow(int iz);
262 : Float_t GetZ2X(int iz);
263 : Float_t GetDZ2X();
264 : Float_t GetDZ2XI();
265 : void FindVoxel(float x, float y2x, float z2x, int &ix,int &ip, int &iz);
266 : void FindVoxel(float x, float y2x, float z2x, UChar_t &ix,UChar_t &ip, UChar_t &iz);
267 : void GetVoxelCoordinates(int isec, int ix, int ip, int iz,float &x, float &p, float &z);
268 : Double_t GetKernelWeight(double *u2vec, int np) const;
269 :
270 : Long64_t GetBin2Fill(const UChar_t binVox[kVoxDim], UShort_t bVal) const;
271 : UShort_t GetVoxGBin(const UChar_t bvox[kVoxDim]) const;
272 : UShort_t GetVoxGBin(int ix,int ip,int iz) const;
273 : void GBin2Vox(UShort_t gbin, UChar_t bvox[kVoxDim]) const;
274 : //
275 0 : void SetRun(int run) {fRun = run;}
276 : void SetTMinMax(Long64_t tmin=0, Long64_t tmax=9999999999);
277 0 : void SetNXBins(int n=kNPadRows) {fNXBins = n;}
278 0 : void SetNY2XBins(int n=15) {fNY2XBins = n;}
279 0 : void SetNZ2XBins(int n=5) {fNZ2XBins = n;}
280 0 : void SetMaxTracks(int n=4000000) {fMaxTracks = n;}
281 0 : void SetNominalTimeBin(int nsec=2400) {fNominalTimeBin = nsec;}
282 0 : void SetNominalTimeBinPrec(float p=0.3) {fNominalTimeBinPrec = p;}
283 0 : void SetFixAligmentBug(Bool_t v=kTRUE) {fFixAlignmentBug = v;}
284 0 : void SetCacheLearnSize(int n=1) {fLearnSize = n;}
285 0 : void SetCacheInput(Int_t v=100) {fCacheInp = v;}
286 0 : void SetSwitchCache(Bool_t v=kFALSE) {fSwitchCache = v;}
287 0 : void SetApplyZt2Zc(Bool_t v=kTRUE) {fApplyZt2Zc = v;}
288 0 : void SetResidualList(const char* l) {fResidualList = l;}
289 0 : void SetOCDBPath(const char* l) {fOCDBPath = l;}
290 0 : void SetUseErrorInSmoothing(Bool_t v=kTRUE) {fUseErrInSmoothing = v;}
291 0 : void SetNPrimTrackCuts(int n=400) {fNPrimTracksCut = n;}
292 0 : void SetMinTrackToUse(int n=600000) {fMinTracksToUse = n;}
293 0 : void SetMinTracksPerVDBin(int n=3000) {fMinTracksPerVDBin = n;}
294 0 : void SetMinEntriesVoxel(int n=15) {fMinEntriesVoxel = n;}
295 0 : void SetMinNClusters(int n=30) {fMinNCl = n;}
296 0 : void SetNVoisinMA(int n=3) {fNVoisinMA = n;}
297 0 : void SetNVoisinMALong(int n=15) {fNVoisinMALong = n;}
298 0 : void SetMaxDevYHelix(float d=0.3) {fMaxDevYHelix = d;}
299 0 : void SetMaxDevZHelix(float d=0.3) {fMaxDevZHelix = d;}
300 0 : void SetMaxStdDevMA(float v=25.0) {fMaxStdDevMA = v;}
301 0 : void SetMaxRMSLong(float v=0.8) {fMaxRMSLong = v;}
302 0 : void SetMaxRejFrac(float v=0.15) {fMaxRejFrac = v;}
303 0 : void SetFilterOutliers(Bool_t v=kTRUE) {fFilterOutliers = v;}
304 :
305 0 : void SetMaxFitYErr2(float v=1.0) {fMaxFitYErr2 = v;}
306 0 : void SetMaxFitXErr2(float v=9.0) {fMaxFitXErr2 = v;}
307 0 : void SetMaxFitXYCorr(float v=0.95) {fMaxFitXYCorr = v;}
308 0 : void SetLTMCut(float v=0.75) {fLTMCut = v;}
309 :
310 0 : Bool_t GetXBinIgnored(int sect, int bin) const {return fXBinIgnore[sect].TestBitNumber(bin);}
311 0 : void SetXBinIgnored(int sect, int bin, Bool_t v=kTRUE) {fXBinIgnore[sect].SetBitNumber(bin,v);}
312 : //
313 0 : Float_t GetMinValidVoxFracDrift() const {return fMinValidVoxFracDrift;}
314 0 : Float_t GetMaxSigY() const {return fMaxSigY;}
315 0 : Float_t GetMaxSigZ() const {return fMaxSigZ;}
316 0 : Int_t GetMaxBadXBinsToCover() const {return fMaxBadXBinsToCover;}
317 0 : Int_t GetMinGoodXBinsToCover() const {return fMinGoodXBinsToCover;}
318 0 : Int_t GetMaxBadRowsPerSector() const {return fMaxBadRowsPerSector;}
319 :
320 0 : void SetMinValidVoxFracDrift(float v=0.7) {fMinValidVoxFracDrift = v;}
321 0 : void SetMaxSigY(Float_t v=1.1) {fMaxSigY = v;}
322 0 : void SetMaxSigZ(Float_t v=0.7) {fMaxSigZ = v;}
323 0 : void SetMaxBadXBinsToCover(int n=4) {fMaxBadXBinsToCover = n;}
324 0 : void SetMinGoodXBinsToCover(int n=2) {fMinGoodXBinsToCover = n;}
325 0 : void SetMaxBadRowsPerSector(float v=0.5) {fMaxBadRowsPerSector = v;}
326 : //
327 0 : Bool_t GetUseTOFBC() const {return fUseTOFBC;}
328 0 : Float_t GetTOFBCMin() const {return fTOFBCMin;}
329 0 : Float_t GetTOFBCMax() const {return fTOFBCMax;}
330 0 : void SetUseTOFBC(Bool_t v) {fUseTOFBC = v;}
331 0 : void SetTOFBCMin(Float_t v=-25.f) {fTOFBCMin = v;}
332 0 : void SetTOFBCMax(Float_t v=50.f) {fTOFBCMax = v;}
333 : //
334 0 : Float_t GetMaxFitYErr2() const {return fMaxFitYErr2;}
335 0 : Float_t GetMaxFitXErr2() const {return fMaxFitXErr2;}
336 0 : Float_t GetMaxFitXYCorr() const {return fMaxFitXYCorr;}
337 0 : Float_t GetLTMCut() const {return fLTMCut;}
338 :
339 0 : Int_t GetRun() const {return fRun;}
340 0 : Long64_t GetTMin() const {return fTMin;}
341 0 : Long64_t GetTMax() const {return fTMax;}
342 0 : Long64_t GetTMinGRP() const {return fTMinGRP;}
343 0 : Long64_t GetTMaxGRP() const {return fTMaxGRP;}
344 0 : Long64_t GetTMinCTP() const {return fTMinCTP;}
345 0 : Long64_t GetTMaxCTP() const {return fTMaxCTP;}
346 0 : Int_t GetNXBins() const {return fNXBins;}
347 0 : Int_t GetNY2XBins() const {return fNY2XBins;}
348 0 : Int_t GetNZ2XBins() const {return fNZ2XBins;}
349 0 : Int_t GetMaxTracks() const {return fMaxTracks;}
350 0 : Int_t GetNominalTimeBin() const {return fNominalTimeBin;}
351 0 : Float_t GetNominalTimeBinPrec() const {return fNominalTimeBinPrec;}
352 0 : Int_t GetCacheInput() const {return fCacheInp;}
353 0 : Int_t GetCacheLearnSize() const {return fLearnSize;}
354 0 : Int_t GetNPrimTrackCuts() const {return fNPrimTracksCut;}
355 0 : Int_t GetMinEntriesVoxel() const {return fMinEntriesVoxel;}
356 0 : Int_t GetMinNClusters() const {return fMinNCl;}
357 0 : Int_t GetNVoisinMA() const {return fNVoisinMA;}
358 0 : Int_t GetNVoisinMALong() const {return fNVoisinMALong;}
359 0 : Float_t GetMaxDevYHelix() const {return fMaxDevYHelix;}
360 0 : Float_t GetMaxDevZHelix() const {return fMaxDevZHelix;}
361 0 : Float_t GetMaxStdDevMA() const {return fMaxStdDevMA;}
362 0 : Float_t GetMaxRMSLong() const {return fMaxRMSLong;}
363 0 : Float_t GetMaxRejFrac() const {return fMaxRejFrac;}
364 0 : Bool_t GetFilterOutliers() const {return fFilterOutliers;}
365 0 : Int_t GetExternalDetectors() const {return fExtDet;}
366 : void SetExternalDetectors(int det=kUseTRDonly);
367 :
368 : //
369 0 : void SetChebZSlicePerSide(int n=1) {fChebZSlicePerSide = n;}
370 0 : void SetChebPhiSlicePerSector(int n=1) {fChebPhiSlicePerSector = n;}
371 0 : Int_t GetChebZSlicePerSide() const {return fChebZSlicePerSide;}
372 0 : Int_t GetChebPhiSlicePerSector() const {return fChebPhiSlicePerSector;}
373 0 : Int_t* GetNPCheb(int dim) const {return (Int_t*)fNPCheb[dim];}
374 0 : Float_t* GetChebPrecD() const {return (Float_t*)fChebPrecD;}
375 : //
376 0 : Bool_t GetFixAlignmentBug() const {return fFixAlignmentBug;}
377 0 : Bool_t GetSwitchCache() const {return fSwitchCache;}
378 0 : Bool_t GetApplyZt2Zc() const {return fApplyZt2Zc;}
379 0 : Bool_t GetUseErrorInSmoothing() const {return fUseErrInSmoothing;}
380 :
381 0 : const TString& GetOCDBPath() const {return fOCDBPath;}
382 0 : const TString& GetReisdualList() const {return fResidualList;}
383 :
384 0 : const AliTPCChebCorr* GetChebCorrObject() const {return fChebCorr;}
385 :
386 : static AliTPCDcalibRes* Load(const char* fname="alitpcdcalibres.root");
387 : static TTree* LoadTree(const char* fname="voxelResTree.root",Int_t mStyle=7,Int_t mCol=kBlack);
388 0 : static void SetUsedInstance(AliTPCDcalibRes* inst) {fgUsedInstance = inst;}
389 0 : static AliTPCDcalibRes* GetUsedInstance() {return fgUsedInstance;}
390 0 : static float GetTPCRowX(int r) {return kTPCRowX[r];}
391 : protected:
392 : //
393 : Bool_t fInitDone; // init flag
394 : Bool_t fUseErrInSmoothing; // weight kernel by point error
395 : Bool_t fSwitchCache; // reset the cache when the reading mode is changing
396 : Bool_t fFixAlignmentBug; // flag to apply the fix
397 : Bool_t fApplyZt2Zc; // Apply fix for using Z_track instead of Z_cluster in the data
398 : // --------------------------------Chebyshev object creation
399 : Int_t fChebZSlicePerSide; // z partitions per side
400 : Int_t fChebPhiSlicePerSector; // azimuthal partitions per sector
401 : Int_t fNPCheb[kResDim][2]; // cheb. nodes per slice
402 :
403 : Float_t fChebPrecD[kResDim]; // nominal precision per output dimension
404 : AliTPCChebCorr* fChebCorr; // final Chebyshev object
405 : // -------------------------------Task defintion
406 : Int_t fRun; // run number
407 : Int_t fExtDet; // external detectors to use
408 : Long64_t fTMin; // time start for timebin
409 : Long64_t fTMax; // time stop for timebin
410 : Long64_t fTMinGRP; // time start from GRP
411 : Long64_t fTMaxGRP; // time stop from GRP
412 : Long64_t fTMinCTP; // time start from CTP
413 : Long64_t fTMaxCTP; // time stop from CTP
414 : Int_t fDeltaTVD; // vdrift T-bin size
415 : //
416 : Int_t fSigmaTVD; // vdrift smoothing window
417 : Int_t fMaxTracks; // max tracks to accept
418 : Int_t fNominalTimeBin;// nominal time bin
419 : Float_t fNominalTimeBinPrec; // allow deviation of time bin from nominal within this limit
420 : Int_t fCacheInp; // input trees cache in MB
421 : Int_t fLearnSize; // event to learn for the cache
422 : Float_t fBz; // B field
423 : Bool_t fDeleteSectorTrees; // delete residuals trees once statistics tree is done
424 : TString fResidualList; // text file with list of residuals tree
425 : TObjArray* fInputChunks; // list of input files used
426 : TString fOCDBPath; // ocdb path
427 : // ------------------------------Selection/filtering cuts
428 : Int_t fMinTracksPerVDBin; // min number of tracks per VDrift bin
429 : Int_t fMinTracksToUse; // produce warning if n tracks is too low
430 : Int_t fMinEntriesVoxel; // min number of entries per voxel to consider
431 : Int_t fNPrimTracksCut; // of >0, cut on event multiplicity
432 : Float_t fMinNCl; // min number of TPC clusters to consider
433 : Float_t fMaxDevYHelix; // max-min Y deviation of interpolating track from helix
434 : Float_t fMaxDevZHelix; // max-min Z deviation of interpolating track from helix
435 : Float_t fNVoisinMA; // N neighbours for moving average
436 : Float_t fNVoisinMALong; // max RMS of cleaned residuals wrt its fNVoisinMALong moving average
437 : Float_t fMaxStdDevMA; // max cluster N std.dev (Y^2+Z^2) wrt moving av. to accept
438 : Float_t fMaxRMSLong; // max RMS of cleaned residuals wrt its fNVoisinMALong moving average
439 : Float_t fMaxRejFrac; // max outlier clusters tagged to accept the track
440 : Float_t fTOFBCMin; // min dTOF cut in ns if validation requested
441 : Float_t fTOFBCMax; // max dTOF cut in ns if validation requested
442 : Bool_t fUseTOFBC; // require TOF BC validation
443 : Bool_t fFilterOutliers; // reject outliers
444 : Bool_t fFatalOnMissingDrift; // if vdrift is needed by absent, produce fatal
445 : Float_t fMaxFitYErr2; // cut on median fit Y err^2
446 : Float_t fMaxFitXErr2; // cut on median fit X err^2
447 : Float_t fMaxFitXYCorr; // cut on max correlation of X,Y errors in median fit
448 : Float_t fLTMCut; // LTM cut for outliers suppression
449 : //
450 : //-------------------------------- voxel validation
451 : Float_t fMaxSigY; // cut on sigY (after slop removal) MAD estimator
452 : Float_t fMaxSigZ; // cut on sigZ LTM estimator
453 : Float_t fMinValidVoxFracDrift; // smooth/parameterize only Xbins with fraction of valid voxels above the threshold
454 : Int_t fMaxBadXBinsToCover; // do not extrapolate to more than this number of bad Xbins
455 : Int_t fMinGoodXBinsToCover; // requre at least this amount of consecutive good bins to parameterize
456 : Float_t fMaxBadRowsPerSector; // block whole sector once the fraction of bad Xbins exceeds
457 :
458 : // -------------------------------Binning
459 : Int_t fNY2XBins; // y/x bins per sector
460 : Int_t fNZ2XBins; // z/x bins per sector
461 : Int_t fNXBins; // n bins in radial dim.
462 : Int_t fNXYBinsProd; // nx*ny bins
463 : Int_t fNBins[kVoxDim]; // bins
464 : Bool_t fUniformBins[kVoxDim]; // uniform binning? Currently only X may be non-uniform (per pad-row)
465 :
466 : Float_t fDZ2X; // Z2X bin size
467 : Float_t fDX; // X bin size
468 : Float_t fDZ2XI; // inverse Z2X bin size
469 : Float_t fDXI; // inverse X bin size
470 :
471 : Int_t fNGVoxPerSector; // total number of geometrical voxels per sector
472 :
473 : Float_t *fMaxY2X; //[fNXBins] max Y/X at each X bin, account for dead zones
474 : Float_t *fDY2X; //[fNXBins] Y/X bin size at given X bin
475 : Float_t *fDY2XI; //[fNXBins] inverse of Y/X bin size at given X bin
476 :
477 : Long64_t fNBProdSt[kVoxHDim]; // aux arrays for fast bin calculation
478 : Int_t fNBProdSectG[kVoxDim]; // aux info for fast bin index calculation in geom voxel space
479 : Int_t fBins[kVoxDim]; // binning in voxel variables
480 :
481 : // ------------------------------Smoothing
482 : Int_t fKernelType; // kernel type
483 : Int_t fStepKern[kVoxDim]; // N bins to consider with given kernel settings
484 : Float_t fKernelWInv[kVoxDim]; // inverse kernel width in bins
485 : Float_t fKernelScaleEdge[kVoxDim]; // optional scaling factors for kernel width on the edge
486 : Bool_t fSmoothPol2[kVoxDim]; // option for use pol1 or pol2 in each direction (no x-terms)
487 : // result of last kernel minimization: value and dV/dX,dV/dY,dV/dZ for each dim
488 : Double_t fLastSmoothingRes[kResDim*kMaxSmtDim]; //! results of last smoothing
489 : // ------------------------------Selection Stats
490 : Int_t fNEvTot; // total number of events
491 : Int_t fNTrSelTot; // selected tracks
492 : Int_t fNTrSelTotWO; // would be selected w/o outliers rejection
493 : Int_t fNReadCallTot; // read calls from input trees
494 : Long64_t fNBytesReadTot; // total bytes read
495 : Int_t fNTestTracks; // number of control tracks in test stat
496 : Float_t fEstTracksPerEvent; // estimate for used tracks per selected event
497 : Float_t fTestStat[kCtrNbr][kCtrNbr]; // control statistics (fraction to fNTestTracks)
498 : TH1F* fEvRateH; // events per time histo
499 : TH1F* fTracksRate; // accepted tracks per second
500 : TH1F* fTOFBCTestH; // TOF BC (ns) for test fNTestTracks tracks
501 : TProfile2D* fHVDTimeInt; // histo used for time integrated vdrift fit
502 : TH2F* fHVDTimeCorr; // histo used for vdrift time correction fit
503 : Float_t fValidFracXBin[kNSect2][kNPadRows]; // fraction of voxels valid per padrow
504 : Int_t fNSmoothingFailedBins[kNSect2]; // number of failed bins/sector, should be 0 to produce parameterization
505 : TBits fXBinIgnore[kNSect2]; // flag to ignore Xbin
506 : Float_t fLumiCOG; // COG lumi for timebin
507 : TGraph* fLumiCTPGraph; // lumi graph from CTP
508 : // ------------------------------VDrift correction
509 : TVectorD *fVDriftParam;
510 : TGraphErrors *fVDriftGraph;
511 : Float_t fCorrTime; //!
512 : // -----------------------------Results of processing
513 : bres_t *fSectGVoxRes[kNSect2]; //! [fNGVoxPerSector] sectors results for geometric voxel
514 : TTree* fStatTree; //! tree with voxels statistics
515 : TTree* fTmpTree[kNSect2]; //! IO tree per sector
516 : TFile* fTmpFile[kNSect2]; //! file for fTmpTree
517 : THnF* fStatHist[kNSect2]; //! histos for statistics bins
518 : TNDArrayT<float> *fArrNDStat[kNSect2]; //! alias arrays for fast access to fStatHist
519 : //
520 : // ----------------------------data exchange structures for trees and between routines
521 : dts_t fDTS; //! binned residuals
522 : dtc_t fDTC; //! corrected residuals for closure test
523 : dtv_t fDTV; //! vdrift residuals
524 : delta_t fDeltaStr; //! input from delta tree
525 : //
526 : // ---------------------------track data-----------------------------------
527 : int fTimeStamp; //! time stamp
528 : int fNCl; //! number of clusters
529 : float fQ2Pt; //! fitted q2pt
530 : float fTgLam; //! fitted tgLambda
531 : float fArrPhi[kNPadRows]; //! cluster phi
532 : float fArrDY[kNPadRows]; //! cluster residual Y
533 : float fArrDZ[kNPadRows]; //! cluster residual Z
534 : float fArrR[kNPadRows]; //! cluster R
535 : float fArrX[kNPadRows]; //! cluster X in sector coordinates (row)
536 : float fArrYCl[kNPadRows]; //! cluster Y in sector coordinates
537 : float fArrZCl[kNPadRows]; //! cluster Z
538 : float fArrYTr[kNPadRows]; //! ref track Y in sector coordinates
539 : float fArrZTr[kNPadRows]; //! ref tracz Z in sector coordinates
540 : float fArrTgSlp[kNPadRows]; //! track inclination at padrow
541 : int fArrSectID[kNPadRows]; //! cluster sector id
542 : //
543 : static AliTPCDcalibRes* fgUsedInstance; //! interface instance to use for parameterization
544 : //
545 : static const float kMaxResid; // max range of distortions, must be <= than the double32_t range of dst_t
546 : static const float kMaxResidZVD; // max range of distortions in VDrift calib, must be <= than the double32_t range of dtv_t
547 : static const float kMaxTgSlp; // max range of tgSlope, must be <= than the double32_t range of dst_t
548 : static const float kSecDPhi;
549 : static const float kMaxQ2Pt;
550 : static const float kMaxGaussStdDev; // don't evaluate gaussian outside this number of StdDevs
551 : // static const float kMaxTgSlp;
552 : // static const float kMaxResid; // max allowed residual
553 : static const float kMinX; // min X to cover
554 : static const float kMaxX; // max X to cover
555 : static const float kMaxZ2X; // max z/x
556 : static const float kZLim[2]; // endcap positions
557 : static const char* kDriftResFileName;
558 : static const char* kLocalResFileName;
559 : static const char* kClosureTestFileName;
560 : static const char* kStatOut;
561 : static const char* kResOut;
562 : static const char* kDriftFileName;
563 : static const float kDeadZone; // dead zone on sector edges in cm
564 : static const float kInvalidR; // to signal invalid R
565 : static const float kInvalidRes; // to signal invalid residual
566 : static const ULong64_t kMByte;
567 : static const Float_t kZeroK; // zero kernel weight
568 : static const char* kControlBr[];
569 : static const char* kVoxName[];
570 : static const char* kResName[];
571 : static const char* kModeName[];
572 : static const char* kExtDetName[];
573 :
574 : static const Float_t kTPCRowX[]; // X of the pad-row
575 : static const Float_t kTPCRowDX[]; // pitch in X
576 :
577 6 : ClassDef(AliTPCDcalibRes,15);
578 : };
579 :
580 : //________________________________________________________________
581 : inline Int_t AliTPCDcalibRes::GetXBinExact(float x)
582 : {
583 : // convert X to bin ID, following pad row widths
584 0 : if (fUniformBins[kVoxX]) {
585 0 : int ix = (x-kMinX)*fDXI;
586 0 : return (ix<0 || ix>=fNXBins) ? -2 : ix;
587 : }
588 0 : else return GetRowID(x);
589 0 : }
590 :
591 : //________________________________________________________________
592 : inline Float_t AliTPCDcalibRes::GetY2X(int ix, int iy)
593 : {
594 : // get Y2X bin center for ix,iy bin
595 0 : return (0.5f+iy)*fDY2X[ix] - fMaxY2X[ix];
596 : }
597 :
598 : //________________________________________________________________
599 : inline Float_t AliTPCDcalibRes::GetY2XLow(int ix, int iy)
600 : {
601 : // get Y2X bin low edge for ix,iy bin
602 0 : return iy*fDY2X[ix] - fMaxY2X[ix];
603 : }
604 :
605 : //________________________________________________________________
606 : inline Float_t AliTPCDcalibRes::GetDY2X(int ix)
607 : {
608 : // get Y2X bin size value for ix bin
609 0 : return fDY2X[ix];
610 : }
611 :
612 : //________________________________________________________________
613 : inline Float_t AliTPCDcalibRes::GetDY2XI(int ix)
614 : {
615 : // get Y2X inverse bin size for ix bin
616 0 : return fDY2XI[ix];
617 : }
618 :
619 : //________________________________________________________________
620 : inline Float_t AliTPCDcalibRes::GetX(int i)
621 : {
622 : // low edge of i-th X bin
623 0 : return (fUniformBins[kVoxX]) ? kMinX+(0.5+i)*fDX : kTPCRowX[i];
624 : }
625 :
626 : //________________________________________________________________
627 : inline Float_t AliTPCDcalibRes::GetXLow(int i)
628 : {
629 : // low edge of i-th X bin
630 0 : return fUniformBins[kVoxX] ? kMinX+i*fDX : kTPCRowX[i] - 0.5*kTPCRowDX[i];
631 : }
632 :
633 : //________________________________________________________________
634 : inline Float_t AliTPCDcalibRes::GetDX(int i)
635 : {
636 : // width of i-th X bin
637 0 : return fUniformBins[kVoxX] ? fDX : kTPCRowDX[i];
638 : }
639 :
640 : //________________________________________________________________
641 : inline Float_t AliTPCDcalibRes::GetDXI(int i)
642 : {
643 : // inverse width of i-th X bin
644 0 : return (fUniformBins[kVoxX]) ? fDXI : 1.f/kTPCRowDX[i];
645 : }
646 :
647 : //________________________________________________________________
648 : inline Int_t AliTPCDcalibRes::GetY2XBinExact(float y2x, int ix)
649 : {
650 : // get exact y2x bin at given x range
651 0 : float bf = ( y2x + fMaxY2X[ix] ) * GetDY2XI(ix);
652 0 : if (bf<0) return -1;
653 0 : else if (bf>=fNY2XBins) return fNY2XBins;
654 0 : return int(bf);
655 0 : }
656 :
657 : //________________________________________________________________
658 : inline Int_t AliTPCDcalibRes::GetY2XBin(float y2x, int ix)
659 : {
660 : // get closest y2x bin at given x range
661 0 : int bf = ( y2x + fMaxY2X[ix] ) * GetDY2XI(ix);
662 0 : if (bf<0) bf = 0;
663 0 : else if (bf>=fNY2XBins) fNY2XBins-1;
664 0 : return bf;
665 : }
666 :
667 : //________________________________________________________________
668 : inline Int_t AliTPCDcalibRes::GetZ2XBinExact(float z2x)
669 : {
670 : // get exact z2x bin at given x range (z2x is positive for clusters not changing the side)
671 0 : float bz = z2x*GetDZ2XI();
672 0 : if (bz>=fNZ2XBins) return -1;
673 0 : if (bz<0) bz = 0; // to account for clusters which moved to wrong side
674 0 : return int(bz);
675 0 : }
676 :
677 : //________________________________________________________________
678 : inline Int_t AliTPCDcalibRes::GetZ2XBin(float z2x)
679 : {
680 : // get closest z2x bin (z2x is positive for clusters not changing the side)
681 0 : int bz = z2x*GetDZ2XI();
682 0 : if (bz<0) bz = 0; // to account for clusters which moved to wrong side
683 0 : return bz<fNZ2XBins ? bz : fNZ2XBins-1;
684 : }
685 :
686 : //________________________________________________________________
687 : inline Float_t AliTPCDcalibRes::GetZ2X(int iz)
688 : {
689 : // get Z2X bin center for iz, !! always positive
690 0 : return (0.5f+iz)*GetDZ2X();
691 : }
692 :
693 : //________________________________________________________________
694 : inline Float_t AliTPCDcalibRes::GetZ2XLow(int iz)
695 : {
696 : // get Z2X bin low edge for iz !! bin positive
697 0 : return iz*GetDZ2X();
698 : }
699 :
700 : //________________________________________________________________
701 : inline Float_t AliTPCDcalibRes::GetDZ2X()
702 : {
703 : // get Z2X bin size value
704 0 : return fDZ2X;
705 : }
706 :
707 : //________________________________________________________________
708 : inline Float_t AliTPCDcalibRes::GetDZ2XI()
709 : {
710 : // get Z2X inverse bin size
711 0 : return fDZ2XI;
712 : }
713 :
714 : //_____________________________________
715 : inline void AliTPCDcalibRes::FindVoxel(float x, float y2x, float z2x, int &ix, int &ip, int &iz)
716 : {
717 : // calculate voxel center sector coordinates (wrt sector)
718 0 : ix = GetXBin(x);
719 0 : ip = GetY2XBin(y2x,ix);
720 0 : iz = GetZ2XBin(z2x);
721 : //
722 0 : }
723 :
724 : //_____________________________________
725 : inline void AliTPCDcalibRes::FindVoxel(float x, float y2x, float z2x, UChar_t &ix, UChar_t &ip, UChar_t &iz)
726 : {
727 : // calculate voxel center sector coordinates (wrt sector)
728 0 : ix = GetXBin(x);
729 0 : ip = GetY2XBin(y2x,ix);
730 0 : iz = GetZ2XBin(z2x);
731 : //
732 0 : }
733 :
734 : //_____________________________________
735 : inline void AliTPCDcalibRes::GetVoxelCoordinates(int isec, int ix, int ip, int iz, float &x, float &p, float &z)
736 : {
737 : // calculate voxel center sector coordinates (wrt sector)
738 0 : x = GetX(ix);
739 0 : p = GetY2X(ix,ip);
740 0 : z = GetZ2X(iz);
741 0 : if (isec>=kNSect) z = -z;
742 0 : }
743 :
744 :
745 : //_____________________________________________________
746 : inline Long64_t AliTPCDcalibRes::GetBin2Fill(const UChar_t binVox[kVoxDim], UShort_t bVal) const
747 : {
748 : // TH4 bin calculation, bval is the last dimention binID
749 0 : ULong64_t binToFill = bVal+1; // 0 bin is undeflow
750 0 : for (int id=kVoxDim;id--;) binToFill += fNBProdSt[id]*(1+binVox[id]);
751 0 : return binToFill;
752 : }
753 :
754 : //_____________________________________________________
755 : inline UShort_t AliTPCDcalibRes::GetVoxGBin(const UChar_t bvox[kVoxDim]) const
756 : {
757 : // index of geometrix voxel
758 0 : int binToFill = bvox[kVoxDim-1];
759 0 : for (int id=kVoxDim-1;id--;) binToFill += fNBProdSectG[id]*bvox[id];
760 0 : return binToFill;
761 : }
762 :
763 : //_____________________________________________________
764 : inline UShort_t AliTPCDcalibRes::GetVoxGBin(int ix,int ip,int iz) const
765 : {
766 : // index of geometrix voxel
767 0 : UChar_t bvox[kVoxDim];
768 0 : bvox[kVoxX] = ix; bvox[kVoxF] = ip; bvox[kVoxZ] = iz;
769 0 : return GetVoxGBin(bvox);
770 0 : }
771 :
772 : //_____________________________________________________
773 : inline void AliTPCDcalibRes::GBin2Vox(UShort_t gbin, UChar_t bvox[kVoxDim]) const
774 : {
775 : // index to geometrix voxel
776 0 : bvox[kVoxDim-1] = gbin%fNBins[kVoxDim-1];
777 0 : for (int id=kVoxDim-1;id--;) bvox[id] = (gbin/fNBProdSectG[id])%fNBins[id];
778 0 : }
779 :
780 : //_____________________________________________
781 : inline void AliTPCDcalibRes::SetTMinMax(Long64_t tmin, Long64_t tmax) {
782 : // set min max time
783 0 : fTMin = tmin;
784 0 : fTMax = tmax;
785 0 : if (fTMin>=fTMax) {
786 0 : fTMin = 0;
787 0 : fTMax = 9999999999;
788 0 : }
789 0 : }
790 :
791 : #endif
|
