Line data Source code
1 : #include "AliTPCDcalibRes.h"
2 : #include "AliCDBPath.h"
3 : #include "AliCDBEntry.h"
4 : #include "AliCDBStorage.h"
5 : #include "AliGRPObject.h"
6 : #include "AliTriggerRunScalers.h"
7 : #include "AliTriggerScalersRecord.h"
8 : #include "AliDAQ.h"
9 : #include "AliTPCcalibDB.h"
10 : #include "AliTPCParam.h"
11 : #include "AliLumiTools.h"
12 : #include <TKey.h>
13 : #include <TF2.h>
14 :
15 : using std::swap;
16 :
17 : // this must be standalone f-n, since the signature is important for Chebyshev training
18 : void trainCorr(int row, float* tzLoc, float* corrLoc);
19 :
20 : // make sure these branches are always connected in InitDeltaFile
21 : const char* AliTPCDcalibRes::kModeName[kNCollectModes] = {"VDrift calibration","Distortions extraction","Closure test"};
22 : const char* AliTPCDcalibRes::kExtDetName[kNExtDetComb] = {"TRDonly","TOFonly","ITSonly","TRD|TOF"};
23 : const char* AliTPCDcalibRes::kControlBr[kCtrNbr] = {"itsOK","trdOK","tofOK","tofBC","nPrimTracks"};
24 : const char* AliTPCDcalibRes::kVoxName[AliTPCDcalibRes::kVoxHDim] = {"z2x","y2x","x","N"};
25 : const char* AliTPCDcalibRes::kResName[AliTPCDcalibRes::kResDim] = {"dX","dY","dZ","Disp"};
26 : const float AliTPCDcalibRes::kMaxResid=20.0f;
27 : const float AliTPCDcalibRes::kMaxResidZVD=40.0f;
28 : const float AliTPCDcalibRes::kMaxTgSlp=2.0;
29 :
30 6 : const float AliTPCDcalibRes::kSecDPhi = 20.f*TMath::DegToRad();
31 : const float AliTPCDcalibRes::kMaxQ2Pt = 3.0f;
32 : const float AliTPCDcalibRes::kMinX = 85.0f;
33 : const float AliTPCDcalibRes::kMaxX = 246.0f;
34 : const float AliTPCDcalibRes::kMaxZ2X = 1.0f;
35 : const float AliTPCDcalibRes::kZLim[2] = {2.49725e+02,2.49698e+02};
36 : const char* AliTPCDcalibRes::kDriftResFileName = "tmpDriftTree";
37 : const char* AliTPCDcalibRes::kLocalResFileName = "tmpDeltaSect";
38 : const char* AliTPCDcalibRes::kClosureTestFileName = "closureTestSect";
39 : const char* AliTPCDcalibRes::kStatOut = "voxelStat";
40 : const char* AliTPCDcalibRes::kResOut = "voxelRes";
41 : const char* AliTPCDcalibRes::kDriftFileName= "fitDrift";
42 : const float AliTPCDcalibRes::kDeadZone = 1.5f;
43 : const float AliTPCDcalibRes::kZeroK = 1e-6;
44 : const float AliTPCDcalibRes::kInvalidR = 10.f;
45 : const float AliTPCDcalibRes::kInvalidRes = -900.0f;
46 : const float AliTPCDcalibRes::kMaxGaussStdDev = 5.0f;
47 : const ULong64_t AliTPCDcalibRes::kMByte = 1024LL*1024LL;
48 :
49 : const Float_t AliTPCDcalibRes::kTPCRowX[AliTPCDcalibRes::kNPadRows] = { // pad-row center X
50 : 85.225, 85.975, 86.725, 87.475, 88.225, 88.975, 89.725, 90.475, 91.225, 91.975, 92.725, 93.475, 94.225, 94.975, 95.725,
51 : 96.475, 97.225, 97.975, 98.725, 99.475,100.225,100.975,101.725,102.475,103.225,103.975,104.725,105.475,106.225,106.975,
52 : 107.725,108.475,109.225,109.975,110.725,111.475,112.225,112.975,113.725,114.475,115.225,115.975,116.725,117.475,118.225,
53 : 118.975,119.725,120.475,121.225,121.975,122.725,123.475,124.225,124.975,125.725,126.475,127.225,127.975,128.725,129.475,
54 : 130.225,130.975,131.725,135.100,136.100,137.100,138.100,139.100,140.100,141.100,142.100,143.100,144.100,145.100,146.100,
55 : 147.100,148.100,149.100,150.100,151.100,152.100,153.100,154.100,155.100,156.100,157.100,158.100,159.100,160.100,161.100,
56 : 162.100,163.100,164.100,165.100,166.100,167.100,168.100,169.100,170.100,171.100,172.100,173.100,174.100,175.100,176.100,
57 : 177.100,178.100,179.100,180.100,181.100,182.100,183.100,184.100,185.100,186.100,187.100,188.100,189.100,190.100,191.100,
58 : 192.100,193.100,194.100,195.100,196.100,197.100,198.100,199.350,200.850,202.350,203.850,205.350,206.850,208.350,209.850,
59 : 211.350,212.850,214.350,215.850,217.350,218.850,220.350,221.850,223.350,224.850,226.350,227.850,229.350,230.850,232.350,
60 : 233.850,235.350,236.850,238.350,239.850,241.350,242.850,244.350,245.850
61 : };
62 : const Float_t AliTPCDcalibRes::kTPCRowDX[AliTPCDcalibRes::kNPadRows] = { // pad-row pitch in X
63 : 0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,
64 : 0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,
65 : 0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,0.750,
66 : 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,
67 : 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,
68 : 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,
69 : 1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,
70 : 1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.500
71 : };
72 :
73 :
74 : AliTPCDcalibRes* AliTPCDcalibRes::fgUsedInstance = 0;
75 :
76 :
77 6 : ClassImp(AliTPCDcalibRes)
78 :
79 : //________________________________________
80 0 : AliTPCDcalibRes::AliTPCDcalibRes(int run,Long64_t tmin,Long64_t tmax,const char* resList) :
81 0 : fInitDone(kFALSE)
82 0 : ,fUseErrInSmoothing(kTRUE)
83 0 : ,fSwitchCache(kFALSE)
84 0 : ,fFixAlignmentBug(kTRUE)
85 0 : ,fApplyZt2Zc(kTRUE)
86 :
87 0 : ,fChebZSlicePerSide(1)
88 0 : ,fChebPhiSlicePerSector(1)
89 0 : ,fChebCorr(0)
90 :
91 0 : ,fRun(run)
92 0 : ,fExtDet(kUseTRDonly)
93 0 : ,fTMin(tmin)
94 0 : ,fTMax(tmax)
95 0 : ,fTMinGRP(0)
96 0 : ,fTMaxGRP(0)
97 0 : ,fTMinCTP(0)
98 0 : ,fTMaxCTP(0)
99 0 : ,fDeltaTVD(120)
100 0 : ,fSigmaTVD(600)
101 0 : ,fMaxTracks(4000000)
102 0 : ,fNominalTimeBin(2400)
103 0 : ,fNominalTimeBinPrec(0.3)
104 0 : ,fCacheInp(100)
105 0 : ,fLearnSize(1)
106 0 : ,fBz(0)
107 0 : ,fDeleteSectorTrees(kFALSE) // set to true for production
108 0 : ,fResidualList(resList)
109 0 : ,fInputChunks(0)
110 0 : ,fOCDBPath()
111 :
112 0 : ,fMinTracksToUse(600000)
113 0 : ,fMinTracksPerVDBin(3000)
114 0 : ,fMinEntriesVoxel(15)
115 0 : ,fNPrimTracksCut(400)
116 0 : ,fMinNCl(30)
117 0 : ,fMaxDevYHelix(0.3)
118 0 : ,fMaxDevZHelix(0.3) // !!! VDrift calib. screas up the Z fit, 0.3 w/o vdrift
119 0 : ,fNVoisinMA(3)
120 0 : ,fNVoisinMALong(15)
121 0 : ,fMaxStdDevMA(25.0)
122 0 : ,fMaxRMSLong(0.8)
123 0 : ,fMaxRejFrac(0.15)
124 0 : ,fTOFBCMin(-25.0)
125 0 : ,fTOFBCMax(50.0)
126 0 : ,fUseTOFBC(kFALSE)
127 0 : ,fFilterOutliers(kTRUE)
128 0 : ,fFatalOnMissingDrift(kTRUE)
129 0 : ,fMaxFitYErr2(1.0)
130 0 : ,fMaxFitXErr2(9.)
131 0 : ,fMaxFitXYCorr(0.95)
132 0 : ,fLTMCut(0.75)
133 : //
134 0 : ,fMaxSigY(1.1)
135 0 : ,fMaxSigZ(0.7)
136 0 : ,fMinValidVoxFracDrift(0.7)
137 0 : ,fMaxBadXBinsToCover(4)
138 0 : ,fMinGoodXBinsToCover(3)
139 0 : ,fMaxBadRowsPerSector(0.4)
140 : //
141 0 : ,fNY2XBins(15)
142 0 : ,fNZ2XBins(5)
143 0 : ,fNXBins(-1)
144 0 : ,fNXYBinsProd(0)
145 0 : ,fDZ2X(0)
146 0 : ,fDX(0)
147 0 : ,fDZ2XI(0)
148 0 : ,fDXI(0)
149 0 : ,fNGVoxPerSector(0)
150 : //
151 0 : ,fMaxY2X(0)
152 0 : ,fDY2X(0)
153 0 : ,fDY2XI(0)
154 :
155 0 : ,fKernelType(kGaussianKernel)
156 :
157 0 : ,fNEvTot(0)
158 0 : ,fNTrSelTot(0)
159 0 : ,fNTrSelTotWO(0)
160 0 : ,fNReadCallTot(0)
161 0 : ,fNBytesReadTot(0)
162 0 : ,fNTestTracks(0)
163 0 : ,fEstTracksPerEvent(0)
164 0 : ,fEvRateH(0)
165 0 : ,fTracksRate(0)
166 0 : ,fTOFBCTestH(0)
167 0 : ,fHVDTimeInt(0)
168 0 : ,fHVDTimeCorr(0)
169 0 : ,fLumiCOG(0)
170 0 : ,fLumiCTPGraph(0)
171 0 : ,fVDriftParam(0)
172 0 : ,fVDriftGraph(0)
173 0 : ,fCorrTime(0)
174 :
175 0 : ,fStatTree(0)
176 :
177 0 : ,fDTS()
178 0 : ,fDTC()
179 0 : ,fDTV()
180 0 : ,fDeltaStr()
181 :
182 0 : ,fTimeStamp(0)
183 0 : ,fNCl(0)
184 0 : ,fQ2Pt(0)
185 0 : ,fTgLam(0)
186 0 : {
187 0 : SetTMinMax(tmin,tmax);
188 0 : for (int i=0;i<kResDim;i++) {
189 0 : for (int j=0;j<2;j++) fNPCheb[i][j] = -1; //determine from job binning// 15;
190 0 : fChebPrecD[i] = 100e-4;
191 : }
192 :
193 0 : memset(fLastSmoothingRes,0,kResDim*4*sizeof(double));
194 0 : memset(fTestStat,0,kCtrNbr*kCtrNbr*sizeof(float));
195 0 : for (int i=kVoxDim;i--;) {
196 0 : fUniformBins[i] = kTRUE;
197 0 : fKernelScaleEdge[i] = 1.0f;
198 0 : fStepKern[i] = 1;
199 0 : fKernelWInv[i] = 0; // calculated later
200 0 : fSmoothPol2[i] = kFALSE;
201 : }
202 0 : fSmoothPol2[kVoxX] = kTRUE;
203 0 : fSmoothPol2[kVoxF] = kTRUE;
204 : //
205 0 : for (int i=kVoxHDim;i--;) fNBProdSt[i] = 0;
206 0 : for (int i=kVoxDim;i--;) fNBProdSectG[i] = 0;
207 : //
208 0 : for (int i=0;i<kNSect2;i++) {
209 0 : fSectGVoxRes[i] = 0;
210 0 : fTmpTree[i] = 0;
211 0 : fStatHist[i] = 0;
212 0 : fArrNDStat[i] = 0;
213 0 : fTmpFile[i] = 0;
214 0 : memset(fValidFracXBin[i],0,kNPadRows*sizeof(float));
215 0 : fNSmoothingFailedBins[i] = 0;
216 : }
217 0 : SetKernelType();
218 0 : }
219 :
220 : //________________________________________
221 : AliTPCDcalibRes::~AliTPCDcalibRes()
222 0 : {
223 : // d-tor
224 0 : delete fChebCorr;
225 0 : delete[] fMaxY2X;
226 0 : delete[] fDY2X;
227 0 : delete[] fDY2XI;
228 0 : delete fVDriftParam;
229 0 : delete fVDriftGraph;
230 0 : for (int i=0;i<kNSect2;i++) {
231 0 : delete fSectGVoxRes[i];
232 0 : delete fStatHist[i];
233 : }
234 0 : delete fTracksRate;
235 0 : delete fTOFBCTestH;
236 0 : delete fHVDTimeInt;
237 0 : delete fHVDTimeCorr;
238 0 : delete fInputChunks;
239 0 : }
240 :
241 : //________________________________________
242 : void AliTPCDcalibRes::SetExternalDetectors(int det)
243 : {
244 : // set external detectos choice
245 0 : if (det<0 || det>=kNExtDetComb) {
246 0 : AliErrorF("Invalid external detector %d, allowed range 0:%d, see header file enum{...kNExtDetComb}",det,kNExtDetComb-1);
247 0 : return;
248 : }
249 0 : fExtDet = det;
250 0 : }
251 :
252 : //________________________________________
253 : void AliTPCDcalibRes::CalibrateVDrift()
254 : {
255 : // run VDrift calibration from residual trees.
256 : // Parameters time fDeltaTVD (binning) and fSigmaTVD (smoothing) can be changed from their
257 : // default values by their setters
258 0 : TStopwatch sw;
259 : const float kSafeMargin = 0.2f;
260 : const float kTofEffLow = 0.25f; // if fraction of TOFbc validated to otherwise useful tracks is below this, abandon TOF
261 0 : if (!fInitDone) Init();
262 : // check available statistics
263 0 : if (!EstimateChunkStatistics() || !CollectDataStatistics()) AliFatal("Cannot process further");
264 : //
265 0 : int nChunks = fInputChunks->GetEntriesFast();
266 0 : float ev2Chunk = fNEvTot/nChunks;
267 0 : float fracMult = fTestStat[kCtrNtr][kCtrNtr];
268 0 : Bool_t useTOFBC = fUseTOFBC;
269 : float statEstBCon=0,statEstBCoff=0;
270 : // estimate with TOF BC selection
271 0 : if (fExtDet==kUseTRDonly) statEstBCon = fTestStat[kCtrBC0][kCtrTRD];
272 0 : else if (fExtDet==kUseTOFonly) statEstBCon = fTestStat[kCtrBC0][kCtrTOF];
273 0 : else if (fExtDet==kUseITSonly) statEstBCon = fTestStat[kCtrBC0][kCtrITS];
274 0 : else if (fExtDet==kUseTRDorTOF) statEstBCon = fTestStat[kCtrBC0][kCtrTOF]; // contribution of TRD is negligable
275 : //
276 0 : if (fExtDet==kUseTRDonly) statEstBCoff = fTestStat[kCtrTRD][kCtrTRD];
277 0 : else if (fExtDet==kUseTOFonly) statEstBCoff = fTestStat[kCtrTOF][kCtrTOF];
278 0 : else if (fExtDet==kUseITSonly) statEstBCoff = fTestStat[kCtrITS][kCtrITS];
279 0 : else if (fExtDet==kUseTRDorTOF) statEstBCoff = fTestStat[kCtrTRD][kCtrTRD]+
280 0 : fTestStat[kCtrTOF][kCtrTOF]-fTestStat[kCtrTRD][kCtrTOF];
281 : //
282 0 : Long64_t tmn = fTMin;
283 0 : Long64_t tmx = fTMax;
284 0 : double duration = tmx - tmn;
285 : //
286 0 : statEstBCon *= fTestStat[kCtrNtr][kCtrNtr]; // loss due to the mult selection
287 0 : statEstBCon *= fNTestTracks*nChunks*(1.-kSafeMargin); // safety margin for losses (outliers etc)
288 : //
289 0 : statEstBCoff *= fTestStat[kCtrNtr][kCtrNtr]; // loss due to the mult selection
290 0 : statEstBCoff *= fNTestTracks*nChunks*(1.-kSafeMargin); // safety margin for losses (outliers etc)
291 : //
292 : float tofOK=0.f,tofTot=0.f;
293 0 : if (fTOFBCTestH && fTOFBCMin<fTOFBCMax) {
294 0 : tofOK = fTOFBCTestH->Integral(fTOFBCTestH->FindBin(fTOFBCMin),fTOFBCTestH->FindBin(fTOFBCMax));
295 0 : tofTot = fTOFBCTestH->Integral(1,fTOFBCTestH->GetNbinsX());
296 0 : if (tofTot>0) tofOK /= tofTot;
297 : }
298 0 : AliInfoF("Fraction of TOFBC tracks in %.1f:%.1f window: %.2f (%.2f loss assumed)",
299 : fTOFBCMin,fTOFBCMax,tofOK,kSafeMargin);
300 : //
301 0 : AliInfoF("Estimated tracks stat.for %.1f min: with(w/o) TOF BC selection %d (%d)",
302 : duration/60.,int(statEstBCon),int(statEstBCoff));
303 : //
304 0 : if (statEstBCoff<10) AliFatal("No reasonable data found");
305 :
306 0 : float tofEff = statEstBCon/statEstBCoff;
307 0 : if (tofEff<kTofEffLow) {
308 0 : AliWarningF("TOFBC selection eff. %.3f is below the threshold %.3f",tofEff,kTofEffLow);
309 : useTOFBC = kFALSE;
310 0 : }
311 : int ntr2BinBCon=0,ntr2BinBCoff=0,nTimeBins=0;
312 0 : float tbinDuration = fNominalTimeBin;
313 : // define time bins
314 0 : int nWantedTracks = (fMaxTracks+fMinTracksToUse)/2;
315 0 : AliInfoF("Nominal requested Ntracks/bin: %d (min:%d max:%d)",nWantedTracks,fMinTracksToUse,fMaxTracks);
316 0 : nTimeBins = 1+int(duration/fNominalTimeBin);
317 0 : tbinDuration = duration/nTimeBins;
318 0 : if (tbinDuration/fNominalTimeBin<(1.-fNominalTimeBinPrec)) {
319 0 : nTimeBins = TMath::Max(1,nTimeBins-1);
320 0 : tbinDuration = duration/nTimeBins;
321 0 : }
322 0 : AliInfoF("Suggest %d time bins of %d s durations",nTimeBins,int(tbinDuration));
323 0 : ntr2BinBCon = fNominalTimeBin/duration*statEstBCon;
324 0 : ntr2BinBCoff = fNominalTimeBin/duration*statEstBCoff;
325 0 : AliInfoF("Estimated Ntracks per time bin: with(w/o) TOF BC: %d(%d)",ntr2BinBCon,ntr2BinBCoff);
326 : // if ( ntr2BinBCon<(nWantedTracks+fMinTracksToUse)/2 ) useTOFBC = kFALSE; // try to keep stat. high
327 0 : if ( ntr2BinBCon<nWantedTracks ) useTOFBC = kFALSE; // try to keep stat. high
328 : //
329 0 : if (!useTOFBC && fUseTOFBC) {
330 0 : AliWarning("Switching OFF requested TOF BC validation due to statistics");
331 0 : fUseTOFBC = kFALSE;
332 0 : }
333 0 : fEstTracksPerEvent = (fUseTOFBC ? statEstBCon:statEstBCoff)/fNEvTot;
334 : //
335 0 : printf("StatInfo.NTBin\t%d\n",nTimeBins);
336 0 : printf("StatInfo.TBinDuration\t%d\n",int(tbinDuration));
337 0 : printf("StatInfo.TOFBCsel\t%d\n",fUseTOFBC);
338 : // select tracks matching to time window and write compact local trees
339 : //
340 0 : CollectData(kVDriftCalibMode);
341 : //
342 0 : FitDrift();
343 0 : sw.Stop();
344 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
345 : //
346 0 : }
347 :
348 : //________________________________________
349 : void AliTPCDcalibRes::FitDrift()
350 : {
351 : // fit v.drift params
352 : const int kUsePoints0 = 1000000; // nominal number of points for 1st estimate
353 : const int kUsePoints1 = 10000000; // nominal number of points for time integrated estimate
354 : const int kNXBins = 200, kNYBins = 100, kNTCorrBins=100;
355 : const int kMinEntriesBin = 100;
356 : const float kMaxTCorr = 0.02; // max time correction for vdrift
357 : const float kDiscardDriftEdge = 5.; // discard min and max drift values with this margin
358 0 : TStopwatch sw; sw.Start();
359 0 : AliSysInfo::AddStamp("FitDrift",0,0,0,0);
360 : //
361 0 : if (!fInitDone) Init();
362 : //
363 0 : TString zresFileName = Form("%s.root",kDriftResFileName);
364 0 : TFile* zresFile = TFile::Open(zresFileName.Data());
365 0 : if (!zresFile) AliFatalF("file %s not found",zresFileName.Data());
366 0 : TString treeName = Form("resdrift");
367 0 : TTree *zresTree = (TTree*) zresFile->Get(treeName.Data());
368 0 : if (!zresTree) AliFatalF("tree %s is not found in file %s",treeName.Data(),zresFileName.Data());
369 : //
370 0 : dtv_t *dtvP = &fDTV;
371 0 : zresTree->SetBranchAddress("dtv",&dtvP);
372 0 : int entries = zresTree->GetEntries();
373 0 : if (!entries) AliFatalF("tree %s in %s has no entries",treeName.Data(),zresFileName.Data());
374 : //
375 0 : float prescale = kUsePoints0/float(entries);
376 0 : int npUse = entries*prescale;
377 : //
378 0 : float *xArr = new Float_t[npUse];
379 0 : float *zArr = new Float_t[npUse];
380 : int nAcc = 0;
381 0 : float dmaxCut = (kZLim[0]+kZLim[1])/2. - kDiscardDriftEdge;
382 0 : float res[2],err[3];
383 0 : for (int i=0;i<entries;i++) {
384 0 : if (gRandom->Rndm()>prescale) continue;
385 0 : zresTree->GetEntry(i);
386 0 : if (/*fDTV.drift<kDiscardDriftEdge ||*/ fDTV.drift>dmaxCut) continue;
387 0 : xArr[nAcc] = fDTV.drift;
388 0 : zArr[nAcc] = fDTV.side>0 ? fDTV.dz : -fDTV.dz;
389 0 : nAcc++;
390 0 : }
391 0 : AliInfoF("Will use %d points out of %d for outliers rejection estimate",nAcc,entries);
392 : //
393 0 : float sigMAD = AliTPCDcalibRes::FitPoly1Robust(nAcc,xArr,zArr,res,err,0.95);
394 0 : if (sigMAD<0) AliFatal("Unbinned fit failed");
395 0 : AliInfoF("Will clean outliers outside of |dz*side-(%.3e+drift*(%.3e))|<3*%.3f band",res[0],res[1],sigMAD);
396 0 : delete[] xArr;
397 0 : delete[] zArr;
398 : //
399 0 : const float outlCut = sigMAD*3.0f;
400 0 : delete fHVDTimeInt;
401 0 : fHVDTimeInt = new TProfile2D("driftTInt","",kNXBins,-dmaxCut,dmaxCut,kNYBins,-kMaxX,kMaxX);
402 0 : fHVDTimeInt->SetXTitle("side*drift");
403 0 : fHVDTimeInt->SetYTitle("ylab");
404 0 : fHVDTimeInt->SetZTitle("#deltaz");
405 0 : fHVDTimeInt->SetDirectory(0);
406 0 : prescale = kUsePoints1/float(entries);
407 0 : npUse = entries*prescale;
408 0 : AliInfoF("Will use %d points out of %d for time-integrated estimate",npUse,entries);
409 0 : for (int i=0;i<entries;i++) {
410 0 : if (gRandom->Rndm()>prescale) continue;
411 0 : zresTree->GetEntry(i);
412 0 : if (/*fDTV.drift<kDiscardDriftEdge ||*/ fDTV.drift>dmaxCut) continue;
413 0 : float dz = fDTV.dz;
414 0 : float sdz = fDTV.side>0 ? dz : -dz;
415 0 : if (TMath::Abs(sdz - (res[0]+res[1]*fDTV.drift)) > outlCut) continue;
416 0 : float sdrift = fDTV.side>0 ? fDTV.drift : -fDTV.drift;
417 : // float d2z = fDTV.drift/kZLim[fDTV.side<0];
418 : float wgh = 1;//160./fDTV.r; // use weight as inverse^2 of radius to increase the ITS contribution
419 0 : fHVDTimeInt->Fill(sdrift, fDTV.ylab, dz, wgh);
420 : //fHVDTimeInt->Fill(sdrift, fDTV.ylab*sdrift/kZLim[fDTV.side<0], sdz); // old way
421 0 : }
422 0 : float zlim = (kZLim[0]+kZLim[1])/2.f;
423 0 : TF2* ftd = new TF2("ftd",Form("[0]*sign(x)+[1]+([2]*y/%.2f+[3])*x",zlim),-250,250,-250,250);
424 : //TF2* ftd = new TF2("ftd","[0]+[1]*sign(x)+[2]*sign(x)*y + [3]*sign(x)*x",-250,250,-250,250); // old fun
425 0 : ftd->SetParameters(res[0],0.,0.,res[1]); // initial values from unbinned fit
426 0 : AliInfoF("Fitting time-integrated vdrift params by %s",ftd->GetTitle());
427 0 : TFitResultPtr rf = fHVDTimeInt->Fit(ftd,"0S");
428 0 : int ndf = rf->Ndf();
429 0 : float chi2 = rf->Chi2();
430 0 : AliInfoF("Fit chi2: %f per %d DOFs -> %f",chi2,ndf,ndf>0 ? chi2/ndf : -1.f);
431 0 : delete fVDriftParam;
432 0 : fVDriftParam = new TVectorD(4);
433 0 : for (int i=4;i--;) (*fVDriftParam)[i] = ftd->GetParameter(i);
434 : //
435 : // time correction
436 0 : delete fHVDTimeCorr;
437 0 : int ntbins = double(fTMax-fTMin)/fDeltaTVD+1;
438 0 : fHVDTimeCorr = new TH2F("driftTCorr","drift time correction",ntbins,fTMin,fTMax+1,kNTCorrBins,-kMaxTCorr,kMaxTCorr);
439 0 : fHVDTimeCorr->SetDirectory(0);
440 0 : fHVDTimeCorr->SetXTitle("time");
441 0 : fHVDTimeCorr->SetYTitle("delta");
442 : //
443 0 : double *vpar = fVDriftParam->GetMatrixArray();
444 0 : for (int i=0;i<entries;i++) {
445 0 : zresTree->GetEntry(i);
446 0 : if (/*fDTV.drift<kDiscardDriftEdge ||*/ fDTV.drift>dmaxCut) continue;
447 0 : float dz = fDTV.side>0 ? fDTV.dz : -fDTV.dz;
448 0 : float sdrift = fDTV.side>0 ? fDTV.drift : -fDTV.drift;
449 0 : float d2z = fDTV.drift/kZLim[fDTV.side<0];
450 0 : Double_t expected = vpar[0]+vpar[1]*fDTV.side + vpar[2]*fDTV.ylab*d2z + vpar[3]*fDTV.drift;
451 0 : dz -= expected;
452 0 : if (TMath::Abs(dz) > outlCut) continue;
453 : float wgh = 1;//160./fDTV.r; // use weight as inverse^2 of radius to increase the ITS contribution
454 0 : fHVDTimeCorr->Fill(fDTV.t, dz/fDTV.drift, d2z*wgh); // ?? why this weight
455 0 : }
456 : //
457 : // check results
458 : //
459 : // extract values
460 0 : TF1* gs = new TF1("gs","gaus",-kMaxTCorr,kMaxTCorr);
461 0 : TObjArray arrH;
462 0 : arrH.SetOwner(kTRUE);
463 0 : fHVDTimeCorr->FitSlicesY(gs,0,-1,0,"QNR",&arrH);
464 0 : TH1* hEnt = fHVDTimeCorr->ProjectionX("hEnt",1,ntbins);
465 0 : TH1* hmean = (TH1*)arrH[1];
466 0 : double tArr[ntbins],dArr[ntbins],deArr[ntbins];
467 : nAcc = 0;
468 0 : for (int i=0;i<ntbins;i++) {
469 0 : if (hEnt->GetBinContent(i+1)<kMinEntriesBin) continue;
470 0 : deArr[nAcc] = hmean->GetBinError(i+1);
471 0 : if (deArr[nAcc]<1e-16) continue;
472 0 : tArr[nAcc] = hmean->GetBinCenter(i+1);
473 0 : dArr[nAcc] = hmean->GetBinContent(i+1);
474 0 : nAcc++;
475 0 : }
476 : //
477 0 : delete fVDriftGraph;
478 0 : fVDriftGraph = new TGraphErrors(ntbins);
479 0 : double resve[2];
480 : int ngAcc = 0;
481 0 : for (int i=0;i<ntbins;i++) {
482 0 : double tQuery = hEnt->GetBinCenter(i);
483 0 : if (!GetSmooth1D(tQuery,resve,nAcc,tArr,dArr,deArr,fSigmaTVD,kGaussianKernel,kFALSE,kTRUE)) {
484 0 : AliWarningF("Failed to smooth at point %d out of %d (T=%d)",i,ntbins,int(tQuery));
485 0 : continue;
486 : }
487 0 : fVDriftGraph->SetPoint(ngAcc,tQuery,resve[0]);
488 0 : fVDriftGraph->SetPointError(ngAcc,0.5,resve[1]);
489 0 : ngAcc++;
490 0 : }
491 0 : for (int i=ntbins-1;i>=ngAcc; i--) fVDriftGraph->RemovePoint(i); // eliminate empty points
492 : //
493 : // flag if TOFBC was used for VDrift extraction
494 0 : if (fUseTOFBC) fVDriftGraph->SetBit(kVDWithTOFBC);
495 0 : else fVDriftGraph->SetBit(kVDWithoutTOFBC);
496 : //
497 0 : delete hEnt;
498 0 : arrH.Delete();
499 : //
500 0 : SaveFitDrift();
501 : //
502 0 : sw.Stop();
503 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
504 0 : AliSysInfo::AddStamp("FitDrift",1,0,0,0);
505 0 : }
506 :
507 : //________________________________________
508 : void AliTPCDcalibRes::ProcessFromDeltaTrees()
509 : {
510 : // process from residual trees
511 0 : TStopwatch sw;
512 : // select tracks matching to time window and write compact local trees
513 0 : EstimateChunkStatistics();
514 : //
515 0 : CollectData(kDistExtractMode);
516 : //
517 0 : if (fNTrSelTot<fMinTracksToUse) {
518 0 : AliErrorF("Low statistics: number of contributing tracks %d, min.requested %d",fNTrSelTot,fMinTracksToUse);
519 0 : TString stopOnLosStat = gSystem->Getenv("stopOnLowStat");
520 0 : if (!stopOnLosStat.IsNull()) {
521 0 : AliInfo("Stop on low statistics requested: abandoning map creation");
522 0 : exit(1);
523 : }
524 : else {
525 0 : AliInfo("No stop on low statistics requested: starting map creation");
526 : }
527 0 : }
528 0 : ProcessFromLocalBinnedTrees();
529 0 : sw.Stop();
530 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
531 : //
532 0 : }
533 :
534 : //________________________________________
535 : void AliTPCDcalibRes::ProcessFromLocalBinnedTrees()
536 : {
537 : // process starting from local binned trees created by CollectData(kDistExtractMode)
538 0 : TStopwatch sw;
539 0 : sw.Start();
540 :
541 : // do per-sector projections and fits
542 0 : ProcessResiduals();
543 : //
544 : // ProcessDispersions();
545 : //
546 0 : CreateCorrectionObject();
547 : //
548 0 : WriteResTree();
549 : //
550 0 : sw.Stop();
551 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
552 0 : }
553 :
554 : //________________________________________
555 : void AliTPCDcalibRes::ReProcessFromResVoxTree(const char* resTreeFile, Bool_t backup)
556 : {
557 : // reprocess starting from the raw data filled from existing resVox tree
558 0 : TStopwatch sw;
559 0 : sw.Start();
560 0 : if (!LoadDataFromResTree(resTreeFile)) return;
561 0 : ReProcessResiduals();
562 : //
563 0 : if (fChebCorr) delete fChebCorr; fChebCorr = 0;
564 : //
565 0 : CreateCorrectionObject();
566 : //
567 0 : if (backup) {
568 0 : TString inps = resTreeFile;
569 0 : if (inps == GetVoxResFileName()) {
570 0 : TString inpsb = resTreeFile;
571 0 : inpsb.ReplaceAll(".root","_1.root");
572 0 : rename(inps.Data(),inpsb.Data());
573 0 : AliInfoF("Input file %s backed up to %s",inps.Data(),inpsb.Data());
574 0 : }
575 0 : }
576 0 : WriteResTree();
577 : //
578 0 : sw.Stop();
579 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
580 0 : }
581 :
582 : //________________________________________
583 : AliTPCDcalibRes* AliTPCDcalibRes::Load(const char* fname)
584 : {
585 : // load AliTPCDcalibRes object from input file
586 0 : TString fnames = fname;
587 0 : if (!fnames.EndsWith(".root")) {
588 0 : if (!fnames.EndsWith("/")) fnames += "/";
589 0 : fnames += "alitpcdcalibres.root";
590 : }
591 0 : TFile* fl = TFile::Open(fnames.Data());
592 0 : if (!fl) {AliErrorClassF("Failed to open %s",fnames.Data()); return 0;}
593 0 : TList* lstk = fl->GetListOfKeys();
594 0 : TIter next(lstk);
595 : TKey* key = 0;
596 : AliTPCDcalibRes* res = 0;
597 0 : while (key=(TKey*)next()) {
598 0 : TString keyt = key->GetTitle();
599 0 : if (keyt == "AliTPCDcalibRes") {res = (AliTPCDcalibRes*)fl->Get(key->GetName()); break;}
600 0 : }
601 0 : if (!res) AliErrorClassF("Did not find AliTPCDcalibRes object in %s",fnames.Data());
602 : return res;
603 : //
604 0 : }
605 :
606 : //________________________________________
607 : TTree* AliTPCDcalibRes::LoadTree(const char* fname, int markerStyle,int color)
608 : {
609 : // load voxResTee from input file
610 0 : TString fnames = fname;
611 0 : if (!fnames.EndsWith(".root")) {
612 0 : if (!fnames.EndsWith("/")) fnames += "/";
613 0 : fnames += "voxelResTree.root";
614 : }
615 0 : TFile* fl = TFile::Open(fnames.Data());
616 0 : if (!fl) {AliErrorClassF("Failed to open %s",fnames.Data()); return 0;}
617 0 : TList* lstk = fl->GetListOfKeys();
618 0 : TIter next(lstk);
619 : TKey* key = 0;
620 : TTree* res = 0;
621 0 : while (key=(TKey*)next()) {
622 0 : TString keyt = key->GetName();
623 0 : if (keyt == "voxRes") {res = (TTree*)fl->Get(key->GetName()); break;}
624 0 : }
625 0 : if (!res) AliErrorClassF("Did not find voxRes tree in %s",fnames.Data());
626 : else {
627 0 : res->SetMarkerStyle(markerStyle);
628 0 : res->SetMarkerColor(color);
629 0 : res->SetLineColor(color);
630 : }
631 : return res;
632 : //
633 0 : }
634 :
635 : //________________________________________
636 : void AliTPCDcalibRes::Save(const char* name)
637 : {
638 : // save itself
639 0 : TString names = name;
640 0 : if (names.IsNull()) {
641 : // names = Form("%s_run%d_%lld_%lld.root",IsA()->GetName(),fRun,fTMin,fTMax);
642 0 : names = Form("%s.root",IsA()->GetName());
643 0 : names.ToLower();
644 : }
645 0 : TFile* flout = TFile::Open(names.Data(),"recreate");
646 0 : this->Write("",TObject::kOverwrite);
647 0 : flout->Close();
648 0 : delete flout;
649 0 : AliInfoF("Saved itself to %s",names.Data());
650 : //
651 0 : }
652 :
653 : //________________________________________
654 : void AliTPCDcalibRes::SaveFitDrift()
655 : {
656 : // save drift parameters in old fitDrift.root format for old way of loading the drift parameters
657 0 : TTreeSRedirector *pcstream = new TTreeSRedirector(Form("%s.root",kDriftFileName),"recreate");
658 0 : (*pcstream)<<"fitTimeStat"<<
659 0 : "grTRDReg.="<<fVDriftGraph<< // time dependent drift correction TRD - regression estimator
660 0 : "paramRobust.="<<fVDriftParam<< // time independent parameters
661 0 : "time0="<<fTMin<<
662 0 : "time1="<<fTMax<<
663 0 : "run="<<fRun<<
664 : "\n";
665 0 : delete pcstream;
666 0 : }
667 :
668 : //==================================================================================
669 : void AliTPCDcalibRes::Init()
670 : {
671 : // do the initialization once
672 : const int kMaxResBins = 0xff;
673 0 : AliSysInfo::AddStamp("ProjStart",0,0,0,0);
674 0 : if (fInitDone) {AliInfo("Init already done"); return;}
675 0 : if (fRun<1) {
676 0 : int run = TString(gSystem->Getenv("runNumber")).Atoi();
677 0 : if (run<1) AliFatal("Run number is neither set nor provided via runNumber env.var");
678 0 : SetRun(run);
679 0 : }
680 : //
681 0 : AliCDBManager* man = AliCDBManager::Instance();
682 0 : if (fOCDBPath.IsNull()) fOCDBPath = "raw://";
683 0 : if (!man->IsDefaultStorageSet()) man->SetDefaultStorage(fOCDBPath);
684 0 : if (man->GetRun()!=fRun) man->SetRun(fRun);
685 : //
686 : // memorize GRP time
687 0 : AliGRPObject* grp = (AliGRPObject*)man->Get(AliCDBPath("GRP/GRP/Data"))->GetObject();
688 0 : fTMinGRP = grp->GetTimeStart();
689 0 : fTMaxGRP = grp->GetTimeEnd();
690 : //
691 0 : if (fUseTOFBC) { // check if TOF is present
692 0 : Int_t activeDetectors = grp->GetDetectorMask();
693 0 : if (!(activeDetectors&AliDAQ::DetectorPattern("TOF"))) {
694 0 : AliWarning("Disabling TOF BC validation since TOF is not in the GRP");
695 0 : fUseTOFBC = kFALSE;
696 0 : }
697 0 : if (fTOFBCMin>=fTOFBCMax) {
698 0 : AliWarningF("Disabling TOF BC validation: inconsistent cuts %.3f:%.3f",fTOFBCMin,fTOFBCMax);
699 0 : fUseTOFBC = kFALSE;
700 0 : }
701 0 : }
702 : // get real data time from CTP scalers
703 0 : AliTriggerRunScalers* scalers = (AliTriggerRunScalers*)man->Get(AliCDBPath("GRP/CTP/Scalers"))->GetObject();
704 0 : Int_t nEntries = scalers->GetScalersRecords()->GetEntriesFast();
705 0 : fTMinCTP = scalers->GetScalersRecord(0 )->GetTimeStamp()->GetSeconds();
706 0 : fTMaxCTP = scalers->GetScalersRecord(nEntries-1)->GetTimeStamp()->GetSeconds();
707 : //
708 0 : if (fTMin<fTMinGRP) fTMin = fTMinGRP;
709 0 : if (fTMax>fTMaxGRP) fTMax = fTMaxGRP;
710 : //
711 0 : AliInfoF("Run time: GRP:%lld:%lld |CTP:%lld:%lld |User:%lld:%lld",fTMinGRP,fTMaxGRP,fTMinCTP,fTMaxCTP,fTMin,fTMax);
712 : //
713 0 : InitFieldGeom();
714 0 : SetName(Form("run%d_%lld_%lld",fRun,fTMin,fTMax));
715 0 : SetTitle(IsA()->GetName());
716 : //
717 : // define boundaries
718 0 : InitBinning();
719 : //
720 0 : LoadVDrift(); //!!!
721 : //
722 : // prepare aux info for stat and residuals histo bin calculation, see doc of TNDArray bin calculation
723 0 : fNBProdSt[kVoxHDim-1] = 1;
724 0 : fNBProdSectG[kVoxDim-1] = 1;
725 0 : for (int i=kVoxHDim-1;i--;) { // +2 to account for under/over-flows
726 0 : fNBProdSt[i] = fNBProdSt[i+1]*(2 + ((i==kVoxDim-1) ? kVoxHDim : fNBins[i+1]));
727 : }
728 0 : for (int i=kVoxDim-1;i--;) {
729 0 : fNBProdSectG[i] = fNBProdSectG[i+1]*fNBins[i+1];
730 : }
731 : //
732 0 : AliSysInfo::AddStamp("Init",0,0,0,0);
733 : //
734 0 : fInitDone = kTRUE;
735 0 : }
736 :
737 : //_____________________________________________________
738 : void AliTPCDcalibRes::CloseTreeFile(TTree* dtree)
739 : {
740 : // close input delta chunk
741 : TFile* fl = 0;
742 0 : TDirectory* dir = dtree->GetDirectory();
743 0 : if (dir) fl = dir->GetFile();
744 0 : delete dtree;
745 0 : if (fl) fl->Close();
746 0 : delete fl;
747 0 : }
748 :
749 : //_____________________________________________________
750 : TTree* AliTPCDcalibRes::InitDeltaFile(const char* name, Bool_t connect, const char* treeName)
751 : {
752 : // init residuals delta file, attach necessary branches
753 : //
754 0 : static delta_t *delta = &fDeltaStr;
755 0 : TString fileNameString(name);
756 0 : if (fileNameString.Contains("alien://") && (!gGrid || (gGrid && !gGrid->IsConnected()))) TGrid::Connect("alien://");
757 0 : TFile* file = TFile::Open(fileNameString.Data());
758 0 : if (!file) {
759 0 : AliErrorF("Cannot open file %s",fileNameString.Data());
760 0 : return 0;
761 : }
762 0 : TTree* tree = (TTree*)file->Get(treeName);
763 0 : if (!tree) {
764 0 : AliErrorF("No tree %s in %s",treeName,fileNameString.Data());
765 0 : delete file; file = 0;
766 0 : return 0;
767 : }
768 : //
769 0 : if (!connect) return tree;
770 : //
771 0 : Bool_t needTRD = fExtDet==kUseTRDorTOF || fExtDet==kUseTRDonly;
772 0 : Bool_t needTOF = fExtDet==kUseTRDorTOF || fExtDet==kUseTOFonly;
773 : //
774 0 : tree->SetCacheLearnEntries(fLearnSize);
775 0 : tree->SetCacheSize(0);
776 0 : tree->SetCacheSize(fCacheInp*kMByte);
777 : //
778 0 : tree->SetBranchStatus("*",kFALSE);
779 0 : tree->SetBranchStatus("timeStamp",kTRUE);
780 0 : tree->SetBranchStatus("itsOK",kTRUE);
781 0 : tree->SetBranchStatus("trdOK",kTRUE);
782 0 : tree->SetBranchStatus("tofOK",kTRUE);
783 0 : tree->SetBranchStatus("vecR.",kTRUE);
784 0 : tree->SetBranchStatus("vecSec.",kTRUE);
785 0 : tree->SetBranchStatus("vecPhi.",kTRUE);
786 0 : tree->SetBranchStatus("vecZ.",kTRUE);
787 0 : tree->SetBranchStatus("track.*",kTRUE);
788 0 : tree->SetBranchStatus("npValid",kTRUE);
789 0 : tree->SetBranchStatus("its0.",kTRUE);
790 0 : tree->SetBranchStatus("its1.",kTRUE);
791 0 : tree->SetBranchStatus("tofBC",kTRUE);
792 0 : tree->SetBranchStatus("nPrimTracks",kTRUE);
793 : //
794 0 : tree->SetBranchAddress("timeStamp",&fDeltaStr.timeStamp);
795 0 : tree->SetBranchAddress("itsOK",&fDeltaStr.itsOK);
796 0 : tree->SetBranchAddress("trdOK",&fDeltaStr.trdOK);
797 0 : tree->SetBranchAddress("tofOK",&fDeltaStr.tofOK);
798 0 : tree->SetBranchAddress("vecR.",&fDeltaStr.vecR);
799 0 : tree->SetBranchAddress("vecSec.",&fDeltaStr.vecSec);
800 0 : tree->SetBranchAddress("vecPhi.",&fDeltaStr.vecPhi);
801 0 : tree->SetBranchAddress("vecZ.",&fDeltaStr.vecZ);
802 0 : tree->SetBranchAddress("track.",&fDeltaStr.param);
803 0 : tree->SetBranchAddress("npValid",&fDeltaStr.npValid);
804 0 : tree->SetBranchAddress("its0.",&fDeltaStr.vecDYITS);
805 0 : tree->SetBranchAddress("its1.",&fDeltaStr.vecDZITS);
806 0 : tree->SetBranchAddress("tofBC",&fDeltaStr.tofBC);
807 0 : tree->SetBranchAddress("nPrimTracks",&fDeltaStr.nPrimTracks);
808 : //
809 0 : if (needTRD) {
810 0 : tree->SetBranchStatus("trd0.",kTRUE);
811 0 : tree->SetBranchStatus("trd1.",kTRUE);
812 0 : tree->SetBranchAddress("trd0.",&fDeltaStr.vecDYTRD);
813 0 : tree->SetBranchAddress("trd1.",&fDeltaStr.vecDZTRD);
814 : }
815 0 : if (needTOF) {
816 0 : tree->SetBranchStatus("tof0.",kTRUE);
817 0 : tree->SetBranchStatus("tof1.",kTRUE);
818 0 : tree->SetBranchAddress("tof0.",&fDeltaStr.vecDYTOF);
819 0 : tree->SetBranchAddress("tof1.",&fDeltaStr.vecDZTOF);
820 : }
821 : //
822 0 : tree->GetEntry(0);
823 : //
824 : return tree;
825 0 : }
826 :
827 : //_____________________________________________________
828 : Bool_t AliTPCDcalibRes::CollectDataStatistics()
829 : {
830 : // scan even trees to get time coverage
831 0 : int nChunks = (!fInputChunks) ? ParseInputList() : fInputChunks->GetEntriesFast();
832 0 : AliInfoF("Collecting event rate from %d chunks",nChunks);
833 0 : if (!fInitDone) Init();
834 0 : int nPrimTrack,timeStamp,tmin=0x7fffffff,tmax=0;
835 0 : int nb = int((fTMaxCTP-fTMinCTP)/10.);
836 0 : fEvRateH = new TH1F("EvRate","EventRate",1+fTMaxGRP-fTMinGRP,fTMinGRP,fTMaxGRP+1);
837 0 : fEvRateH->SetDirectory(0);
838 0 : fEvRateH->GetXaxis()->SetTimeFormat();
839 0 : fNEvTot = 0;
840 0 : for (int ichunk=0;ichunk<nChunks;ichunk++) {
841 0 : TTree *tree = InitDeltaFile(fInputChunks->At(ichunk)->GetName(),kFALSE,"eventInfo");
842 0 : if (!tree) continue;
843 0 : TBranch* brt = tree->GetBranch("timeStamp");
844 0 : TBranch* brn = tree->GetBranch("nPrimTrack");
845 0 : if (!brt || !brn) {
846 0 : AliError("Did not find branch timeStamp or nPromTrack in event info tree");
847 0 : CloseTreeFile(tree);
848 0 : return kFALSE;
849 : }
850 0 : brt->SetAddress(&timeStamp);
851 0 : brn->SetAddress(&nPrimTrack);
852 0 : int nent = tree->GetEntries();
853 0 : for (int i=0;i<nent;i++) {
854 0 : brt->GetEntry(i);
855 0 : brn->GetEntry(i);
856 0 : if (fNPrimTracksCut>0 && nPrimTrack>fNPrimTracksCut) continue;
857 0 : fEvRateH->Fill(timeStamp);
858 0 : if (tmin>timeStamp) tmin = timeStamp;
859 0 : if (tmax<timeStamp) tmax = timeStamp;
860 : }
861 0 : fNEvTot += nent;
862 0 : CloseTreeFile(tree);
863 0 : }
864 : //
865 0 : if (fTMin<tmin) fTMin = tmin;
866 0 : if (fTMax>tmax) fTMax = tmax;
867 0 : AliInfoF("%d selected events in %d chunks covering %d<T<%d",fNEvTot,nChunks,tmin,tmax);
868 0 : printf("StatInfo.minTime\t%lld\n",fTMin);
869 0 : printf("StatInfo.maxTime\t%lld\n",fTMax);
870 0 : return kTRUE;
871 0 : }
872 :
873 : //_____________________________________________________
874 : Int_t AliTPCDcalibRes::ParseInputList()
875 : {
876 : // convert text file to array of file names
877 0 : if (fInputChunks) delete fInputChunks;
878 0 : if (fResidualList.IsNull()) {
879 0 : AliError("Input residuals list is not assigned");
880 0 : return 0;
881 : }
882 0 : TString chunkList = gSystem->GetFromPipe(TString::Format("cat %s",fResidualList.Data()).Data());
883 0 : fInputChunks = chunkList.Tokenize("\n");
884 0 : return fInputChunks->GetEntriesFast();
885 0 : }
886 :
887 : //_____________________________________________________
888 : Bool_t AliTPCDcalibRes::EstimateChunkStatistics()
889 : {
890 : // make rough estimate of statistics with different options
891 0 : AliInfo("Performing rough statistics check");
892 0 : if (!fInitDone) Init();
893 0 : if (!AliGeomManager::GetGeometry()) InitFieldGeom(); // in case started from saved object
894 : //
895 : // pick 1st chunk
896 0 : int nChunks = (!fInputChunks) ? ParseInputList() : fInputChunks->GetEntriesFast();
897 : //
898 : TTree *tree = 0;
899 : int chunk=0; // read 1st accessible chunk
900 0 : while ( !(tree=InitDeltaFile(fInputChunks->At(chunk)->GetName())) && ++chunk<nChunks) {}
901 0 : if (!tree) {
902 0 : AliError("No data chunk was accessible");
903 0 : return kFALSE;
904 : }
905 : //
906 : // make sure these branches are always connected in InitDeltaFile
907 : const char* kNeedBr[]={"itsOK","trdOK","tofOK","tofBC","nPrimTracks"};
908 0 : TBranch* br[kCtrNbr];
909 0 : for (int i=0;i<kCtrNbr;i++) {
910 0 : br[i] = tree->GetBranch(kControlBr[i]);
911 0 : if (!br[i]) AliFatalF("Control branch %s is not in the delta tree",kControlBr[i]);
912 0 : if (!br[i]->GetAddress()) AliFatalF("Control branch %s address is not set",kControlBr[i]);
913 0 : if (!tree->GetBranchStatus(kControlBr[i])) AliFatalF("Control branch %s is not active",kControlBr[i]);
914 : }
915 0 : fNTestTracks = tree->GetEntries();
916 0 : memset(fTestStat,0,kCtrNbr*kCtrNbr*sizeof(float));
917 0 : if (fTOFBCTestH) delete fTOFBCTestH;
918 0 : fTOFBCTestH = new TH1F("TOFBCtest",Form("TOF BC for %d tracks",fNTestTracks),1000,-500.,500.);
919 0 : fTOFBCTestH->SetDirectory(0);
920 0 : Bool_t condOK[kCtrNbr];
921 0 : for (int itr=0;itr<fNTestTracks;itr++) {
922 0 : for (int ib=kCtrNbr;ib--;) br[ib]->GetEntry(itr);
923 0 : condOK[kCtrITS] = fDeltaStr.itsOK;
924 0 : condOK[kCtrTRD] = fDeltaStr.trdOK;
925 0 : condOK[kCtrTOF] = fDeltaStr.tofOK;
926 0 : condOK[kCtrBC0] = fDeltaStr.tofBC>fTOFBCMin && fDeltaStr.tofBC<=fTOFBCMax;
927 0 : condOK[kCtrNtr] = fNPrimTracksCut<0 || fDeltaStr.nPrimTracks<=fNPrimTracksCut;
928 0 : for (int i=kCtrNbr;i--;) for (int j=kCtrNbr;j--;) if (condOK[i]&&condOK[j]) fTestStat[i][j]++;
929 0 : fTOFBCTestH->Fill(fDeltaStr.tofBC);
930 : }
931 0 : if (fNTestTracks) for (int i=kCtrNbr;i--;) for (int j=kCtrNbr;j--;) fTestStat[i][j] /= fNTestTracks;
932 0 : AliInfoF("Accepted statistics wrt %d tracks in chunk id=%d (out of %d)",fNTestTracks,chunk,nChunks);
933 0 : printf(" %11s",""); for (int i=0;i<kCtrNbr;i++) printf(" %11s",kControlBr[i]); printf("\n");
934 0 : for (int i=0;i<kCtrNbr;i++) {
935 0 : printf("*%11s",kControlBr[i]);
936 0 : for (int j=0;j<kCtrNbr;j++) printf(" %11.5f",fTestStat[i][j]);
937 0 : printf("\n");
938 : }
939 : //
940 0 : CloseTreeFile(tree);
941 : return kTRUE;
942 0 : }
943 :
944 : //_____________________________________________________
945 : void AliTPCDcalibRes::CollectData(const int mode)
946 : {
947 : const float kEps = 1e-6;
948 : const float q2ptIniTolerance = 1.5;
949 0 : if (!fInitDone) Init();
950 0 : if (!AliGeomManager::GetGeometry()) InitFieldGeom(); // in case started from saved object
951 : //
952 0 : TStopwatch swTot;
953 0 : swTot.Start();
954 0 : fNTrSelTot = 0;
955 0 : fNTrSelTotWO = 0;
956 0 : fNReadCallTot = 0;
957 0 : fNBytesReadTot = 0;
958 : //
959 0 : AliInfo( "***");
960 0 : AliInfoF("*** Collecting data in mode: %s",kModeName[mode]);
961 0 : AliInfoF("*** Ext.Det.used: %s, TOFBC validation: %s",kExtDetName[fExtDet],fUseTOFBC?"ON":"OFF");
962 0 : AliInfo( "***");
963 : //
964 0 : delete fTracksRate;
965 : // init histo for track rate
966 0 : fTracksRate = new TH1F("TracksRate","TracksRate", 1+fTMax-fTMin, fTMin,fTMax+1);
967 0 : fTracksRate->SetDirectory(0);
968 : //
969 0 : delete fTracksRate;
970 : // init histo for track rate
971 0 : fTracksRate = new TH1F("TracksRate","TracksRate", 1+fTMax-fTMin, fTMin,fTMax+1);
972 0 : fTracksRate->SetDirectory(0);
973 : //
974 : float maxAbsResid = kMaxResid - kEps; // discard residuals exceeding this
975 : Bool_t correctVDrift = kTRUE;
976 0 : if (mode==kVDriftCalibMode) {
977 0 : AliInfo("VDrift calibration mode: drift correction disabled");
978 : maxAbsResid = kMaxResidZVD - kEps; // vdrift calibration may see larger residuals
979 : correctVDrift = kFALSE;
980 0 : }
981 : else {
982 0 : if (!fVDriftGraph || !fVDriftParam) {
983 0 : AliErrorF("We are in mode %d but DriftGraph=%p or VDriftParam=%p missing",
984 : mode,fVDriftGraph,fVDriftParam);
985 0 : if (fFatalOnMissingDrift) AliFatal("Aborting...");
986 : }
987 : // make sure TOF usage is consistent with that in VDrift extraction
988 0 : Bool_t vdtofUsed1 = fVDriftGraph->TestBit(kVDWithTOFBC);
989 0 : Bool_t vdtofUsed0 = fVDriftGraph->TestBit(kVDWithoutTOFBC);
990 0 : if (vdtofUsed1 || vdtofUsed0) { // override to what was used in vdrift extraction
991 0 : AliInfoF("VDrift has TOFBC flag %s, imposing in distortion map extraction",vdtofUsed1 ? "ON":"OFF");
992 0 : if (vdtofUsed1) fUseTOFBC = kTRUE;
993 0 : else fUseTOFBC = kFALSE;
994 : }
995 0 : else AliInfoF("VDrift has TOFBC flags set, using user provided setting %s",fUseTOFBC ? "ON":"OFF");
996 : }
997 0 : AliInfoF("TOFBC usage in window %.1f:%.1f is %s",fTOFBCMin,fTOFBCMax,fUseTOFBC ? "ON":"OFF");
998 0 : CreateLocalResidualsTrees(mode);
999 : //
1000 : // if cheb object is present, do on-the-fly init to attach internal structures
1001 0 : if (fChebCorr) fChebCorr->Init();
1002 : // prepare input tree
1003 0 : int nChunks = (!fInputChunks) ? ParseInputList() : fInputChunks->GetEntriesFast();
1004 : //
1005 0 : AliSysInfo::AddStamp("ProjInit",0,0,0,0);
1006 : //
1007 0 : for (int ichunk=0;ichunk<nChunks;ichunk++) {
1008 : //
1009 : int ntrSelChunkWO=0, ntrSelChunk=0,nReadCallsChunk=0,nBytesReadChunk=0;
1010 : //
1011 0 : TStopwatch swc;
1012 0 : swc.Start();
1013 0 : TString deltaFName = fInputChunks->At(ichunk)->GetName();
1014 0 : TTree *tree = InitDeltaFile(deltaFName.Data());
1015 0 : if (!tree) continue;
1016 : //
1017 0 : TBranch* brTime = tree->GetBranch("timeStamp");
1018 0 : TBranch* brTRDOK = tree->GetBranch("trdOK");
1019 0 : TBranch* brTOFOK = tree->GetBranch("tofOK");
1020 0 : TBranch* brITSOK = tree->GetBranch("itsOK");
1021 : TBranch* brTOFBC = 0;
1022 0 : if (fUseTOFBC) brTOFBC = tree->GetBranch("tofBC");
1023 : //
1024 0 : int nTracks = tree->GetEntries();
1025 0 : AliInfoF("Processing %d tracks of %s",nTracks,deltaFName.Data());
1026 :
1027 0 : float residHelixY[kNPadRows],residHelixZ[kNPadRows];
1028 : //
1029 : // reset the cache when swithching between the timeStamp and Event read modes
1030 : Bool_t lastReadMatched = kFALSE;
1031 0 : for (int itr=0;itr<nTracks;itr++) {
1032 0 : nBytesReadChunk += brTime->GetEntry(itr);
1033 0 : fTimeStamp = fDeltaStr.timeStamp;
1034 0 : if (fTimeStamp<fTMin || fTimeStamp>fTMax) {
1035 0 : if (lastReadMatched && fSwitchCache) { // reset the cache
1036 0 : tree->SetCacheSize(0);
1037 0 : tree->SetCacheSize(fCacheInp*kMByte);
1038 : lastReadMatched = kFALSE;
1039 0 : }
1040 : continue;
1041 : }
1042 0 : if (mode==kVDriftCalibMode && EnoughStatForVDrift(fTimeStamp)) continue; // matching but is it needed?
1043 : //
1044 0 : brITSOK->GetEntry(itr);
1045 0 : brTRDOK->GetEntry(itr);
1046 0 : brTOFOK->GetEntry(itr);
1047 0 : if (!fDeltaStr.itsOK) continue;
1048 0 : if (!fDeltaStr.trdOK && fExtDet==kUseTRDonly) continue;
1049 0 : if (!fDeltaStr.tofOK && fExtDet==kUseTOFonly) continue;
1050 0 : if (!fDeltaStr.tofOK && !fDeltaStr.trdOK && fExtDet!=kUseITSonly) continue;
1051 : //
1052 0 : if (brTOFBC && brTOFBC->GetEntry(itr) && (fDeltaStr.tofBC<fTOFBCMin || fDeltaStr.tofBC>fTOFBCMax)) continue;
1053 : //
1054 0 : if (!lastReadMatched && fSwitchCache) { // reset the cache before switching to event reading mode
1055 0 : tree->SetCacheSize(0);
1056 0 : tree->SetCacheSize(fCacheInp*kMByte);
1057 : }
1058 : lastReadMatched = kTRUE;
1059 0 : nBytesReadChunk += tree->GetEntry(itr);
1060 0 : if (fNPrimTracksCut>0 && fDeltaStr.nPrimTracks>fNPrimTracksCut) continue;
1061 : //
1062 0 : fQ2Pt = fDeltaStr.param->GetParameter()[4];
1063 0 : fTgLam = fDeltaStr.param->GetParameter()[3];
1064 0 : if (TMath::Abs(fQ2Pt)>kMaxQ2Pt*q2ptIniTolerance) continue;
1065 : //
1066 0 : const Float_t *vSec= fDeltaStr.vecSec->GetMatrixArray();
1067 0 : const Float_t *vPhi= fDeltaStr.vecPhi->GetMatrixArray();
1068 0 : const Float_t *vR = fDeltaStr.vecR->GetMatrixArray();
1069 0 : const Float_t *vZ = fDeltaStr.vecZ->GetMatrixArray();
1070 0 : const Float_t *vDYITS = fDeltaStr.vecDYITS->GetMatrixArray();
1071 0 : const Float_t *vDZITS = fDeltaStr.vecDZITS->GetMatrixArray();
1072 : //
1073 : const Float_t *vDY=0,*vDZ = 0;
1074 0 : if (fExtDet==kUseTRDonly || (fExtDet==kUseTRDorTOF && fDeltaStr.trdOK)) {
1075 0 : vDY = fDeltaStr.vecDYTRD->GetMatrixArray();
1076 0 : vDZ = fDeltaStr.vecDZTRD->GetMatrixArray();
1077 0 : }
1078 0 : else if (fExtDet==kUseITSonly) { // ignore other detectos
1079 0 : vDY = fDeltaStr.vecDYITS->GetMatrixArray();
1080 0 : vDZ = fDeltaStr.vecDZITS->GetMatrixArray();
1081 0 : }
1082 : else { // only TOF
1083 0 : vDY = fDeltaStr.vecDYTOF->GetMatrixArray();
1084 0 : vDZ = fDeltaStr.vecDZTOF->GetMatrixArray();
1085 : }
1086 : //
1087 0 : fCorrTime = (correctVDrift && fVDriftGraph) ? fVDriftGraph->Eval(fTimeStamp):0; // for VDrift correction
1088 : //
1089 0 : fNCl = 0;
1090 : // 1st iteration: collect data in cluster frame
1091 0 : for (int ip=0;ip<fDeltaStr.npValid;ip++) { // 1st fill selected track data to buffer for eventual outlier rejection
1092 0 : if (vR[ip]<kInvalidR || vDY[ip]<kInvalidRes || vDYITS[ip]<kInvalidRes) continue;
1093 : //
1094 0 : fArrX[fNCl] = -1;
1095 0 : fArrR[fNCl] = vR[ip]; // X (R) is the same for cluster and track
1096 0 : fArrZTr[fNCl] = vZ[ip]; // Z of ITS track was stored!!
1097 0 : fArrDY[fNCl] = vDY[ip]; // this is also the track coordinate in cluster frame
1098 0 : fArrDZ[fNCl] = vDZ[ip];
1099 0 : fArrPhi[fNCl] = vPhi[ip];
1100 0 : int rocID = TMath::Nint(vSec[ip]);
1101 : //
1102 : // !!! fArrZTr corresponds to ITS track Z, we need that of TRD-ITS
1103 0 : fArrZTr[fNCl] += fArrDZ[fNCl] - vDZITS[ip]; // recover ITS-TRD track position from ITS and deltas
1104 :
1105 0 : if (fFixAlignmentBug && !fDeltaStr.param->TestBit(kAlignmentBugFixedBit)) {
1106 0 : FixAlignmentBug(rocID, fQ2Pt, fBz, fArrPhi[fNCl], fArrR[fNCl], fArrZTr[fNCl], fArrDY[fNCl],fArrDZ[fNCl]);
1107 : }
1108 0 : if (fArrPhi[fNCl]<0) fArrPhi[fNCl] += 2.*TMath::Pi();
1109 : //
1110 : // correct for drift velocity calibration if needed
1111 0 : if (correctVDrift) fArrDZ[fNCl] += GetDriftCorrection(fArrZTr[fNCl],fArrR[fNCl],fArrPhi[fNCl],rocID);
1112 : //
1113 0 : fArrSectID[fNCl] = rocID%kNSect2; // 0-36 for sectors from A0 to C17
1114 : //
1115 0 : fNCl++;
1116 0 : }
1117 : // fit track coordinates by helix to get interpolated track q/pt:
1118 : // more precise than the distorted TPC q/pt
1119 0 : if (fNCl<fMinNCl) continue;
1120 : //
1121 0 : ntrSelChunkWO++;
1122 : //
1123 : float q2ptTPC = fQ2Pt;
1124 0 : Bool_t resH = CompareToHelix(residHelixY,residHelixZ);
1125 : //
1126 0 : if (fFilterOutliers && !resH) continue; // too strong deviation to helix, discard track
1127 0 : if (TMath::Abs(fQ2Pt)>kMaxQ2Pt) continue; // now we have more precise estimate of q/pt
1128 : //
1129 : // 2nd iteration: convert everything to sector frame
1130 : // *****************************************************************
1131 : //
1132 : // All these manipulations are needed because the ResidualTree is stored
1133 : // in cluster frame, while we need the sector frame
1134 : //
1135 : // *****************************************************************
1136 0 : int nc0 = fNCl;
1137 0 : fNCl = 0;
1138 0 : for (int ip=0;ip<nc0;ip++) {
1139 0 : int side = ((fArrSectID[ip] /kNSect)&0x1);
1140 0 : float sna = TMath::Sin(fArrPhi[ip]-(0.5f +fArrSectID[ip]%kNSect)*kSecDPhi);
1141 0 : float csa = TMath::Sqrt((1.f-sna)*(1.f+sna));
1142 : //
1143 : // by using propagation in cluster frame in AliTPCcalibAlignInterpolation::Process,
1144 : // the X of the track is evaluated not at the pad-row x=r*csa but at x=r*sca-dy*sna
1145 0 : double xrow = fArrR[ip]*csa;
1146 0 : double dx = fArrDY[ip]*sna;
1147 : double xtr = xrow - dx;
1148 0 : double ycl = fArrR[ip]*sna; // cluster Y in the sector frame
1149 0 : double ytr = ycl + fArrDY[ip]*csa; // track Y in the sector frame at x=xtr is
1150 : //
1151 0 : double ztr = fArrZTr[ip]; // Z of the track at x=xtr
1152 0 : double zcl = ztr - fArrDZ[ip]; // and the Z of the cluster is Ztr-deltaZ
1153 : //
1154 : // Now we need to take the track to real pad-row X
1155 : // use linear extrapolation:
1156 0 : float tgs = fArrTgSlp[ip];
1157 0 : if (TMath::Abs(tgs)>kMaxTgSlp) continue;
1158 0 : ytr += dx*tgs;
1159 0 : double csXtrInv = TMath::Sqrt(1.+tgs*tgs); // (inverse cosine of track angle)
1160 0 : ztr += dx*fTgLam*csXtrInv;
1161 : //
1162 : // assign to arrays and recalculate residuals
1163 0 : fArrX[fNCl] = xrow;
1164 0 : fArrYTr[fNCl] = ytr;
1165 0 : fArrZTr[fNCl] = ztr;
1166 : //
1167 0 : fArrYCl[fNCl] = ycl;
1168 0 : fArrZCl[fNCl] = zcl;
1169 0 : fArrDY[fNCl] = ytr - ycl;
1170 0 : fArrDZ[fNCl] = ztr - zcl;
1171 : //
1172 : // we don't want under/overflows
1173 0 : if (TMath::Abs(fArrDY[fNCl])>maxAbsResid) continue;
1174 0 : if (TMath::Abs(fArrDZ[fNCl])>maxAbsResid) continue;
1175 : //
1176 0 : if (fArrX[fNCl]<kMinX || fArrX[fNCl]>kMaxX) continue;
1177 0 : if (TMath::Abs(fArrZCl[fNCl])>kZLim[side]) continue;;
1178 : //
1179 : // End of manipulations to go to the sector frame
1180 : //
1181 0 : fNCl++;
1182 0 : }
1183 :
1184 0 : if (fFilterOutliers && !ValidateTrack()) continue;
1185 :
1186 0 : ntrSelChunk++;
1187 :
1188 0 : switch(mode) {
1189 0 : case kVDriftCalibMode: FillDriftResidualsTrees(); break;
1190 0 : case kDistExtractMode: FillLocalResidualsTrees(); break;
1191 0 : case kDistClosureTestMode: FillCorrectedResiduals(); break;
1192 0 : default: AliFatalF("Uknown mode %d",mode);
1193 0 : };
1194 0 : } // loop over tracks
1195 : //
1196 0 : swc.Stop();
1197 0 : TFile* chunkFile = tree->GetDirectory()?tree->GetDirectory()->GetFile():0;
1198 0 : nReadCallsChunk = chunkFile ? chunkFile->GetReadCalls():0;
1199 0 : AliInfoF("Chunk%3d: selected %d tracks (%d with outliers) from chunk %d | %.1f MB read in %d read calls",
1200 0 : ichunk,ntrSelChunk,ntrSelChunkWO, ichunk,float(nBytesReadChunk)/kMByte,nReadCallsChunk); swc.Print();
1201 0 : fNTrSelTot += ntrSelChunk;
1202 0 : fNTrSelTotWO += ntrSelChunkWO;
1203 0 : fNReadCallTot += nReadCallsChunk;
1204 0 : fNBytesReadTot += nBytesReadChunk;
1205 : //
1206 0 : CloseTreeFile(tree);
1207 0 : AliSysInfo::AddStamp("ProjTreeLoc", ichunk ,fNTrSelTot,fNTrSelTot,fNReadCallTot );
1208 : //
1209 0 : if (fNTrSelTot > fMaxTracks) {
1210 0 : AliInfo("Max number of tracks exceeded");
1211 0 : break;
1212 : }
1213 0 : if (mode==kVDriftCalibMode && EnoughStatForVDrift()) {
1214 0 : AliInfo("Sifficient statistics for VDrift calibration collected");
1215 0 : break;
1216 : }
1217 : //
1218 0 : } // loop over chunks
1219 : //
1220 : // write/close local trees
1221 0 : CloseLocalResidualsTrees(mode);
1222 : //
1223 0 : AliInfoF("Summary: selected %d tracks (%d with outliers) | %.1f MB read in %d read calls",
1224 : fNTrSelTot,fNTrSelTotWO,float(fNBytesReadTot)/kMByte,fNReadCallTot);
1225 0 : swTot.Print();
1226 :
1227 0 : AliSysInfo::AddStamp("ProjTreeLocSave");
1228 :
1229 0 : if (mode==kDistExtractMode) WriteStatHistos();
1230 : //
1231 0 : }
1232 :
1233 : //________________________________________________
1234 : Bool_t AliTPCDcalibRes::EnoughStatForVDrift(int tstamp, float maxHolesFrac)
1235 : {
1236 : // check if collected stat. is enough for VDrift calibration
1237 : // when called with positive timestamp, check just the occupancy of the bin it belong to
1238 : static Bool_t *binStat = 0;
1239 : static int nBinStat = 0;
1240 0 : if (!binStat) {
1241 0 : nBinStat = (fTMax-fTMin)/fDeltaTVD+1;
1242 0 : binStat = new Bool_t[nBinStat];
1243 0 : memset(binStat,0,nBinStat*sizeof(Bool_t));
1244 0 : }
1245 : //
1246 0 : int binv = (tstamp-fTMin)/fDeltaTVD; // time stamp provided: query for specific time bin
1247 0 : if (tstamp>0) {
1248 0 : if (binv>=0 && binv<nBinStat) return binStat[binv];
1249 0 : return kTRUE; // this should not happen normally
1250 : }
1251 : //
1252 : // query for the whole run
1253 : int nHoles = 0, nChecked = 0;
1254 0 : int nbinsT = fTracksRate->GetNbinsX();
1255 0 : int nbinsE = fEvRateH->GetNbinsX();
1256 0 : int ntrCumul[nbinsT+2], sumNtr=0;
1257 0 : int nevCumul[nbinsE+2], sumNev=0;
1258 0 : ntrCumul[0] = nevCumul[0] = 0;
1259 0 : for (int ib=1;ib<=nbinsT;ib++) ntrCumul[ib] = (sumNtr += fTracksRate->GetBinContent(ib));
1260 0 : for (int ib=1;ib<=nbinsE;ib++) nevCumul[ib] = (sumNev += fEvRateH->GetBinContent(ib));
1261 0 : ntrCumul[nbinsT+1] = sumNtr;
1262 0 : nevCumul[nbinsE+1] = sumNev;
1263 : //
1264 : // mean number of tracks expected per tested bin in full stat
1265 0 : float nexpTracksAv = fEstTracksPerEvent*fNEvTot*fDeltaTVD/(fTMax-fTMin);
1266 : int longestEmpty=0,prevEmpty=0;
1267 0 : for (int ts=fTMin;ts<fTMax-fDeltaTVD;ts+=fDeltaTVD) {
1268 0 : int bminT = fTracksRate->FindBin(ts);
1269 0 : int bmaxT = fTracksRate->FindBin(ts+fDeltaTVD);
1270 0 : int bminE = fEvRateH->FindBin(ts);
1271 0 : int bmaxE = fEvRateH->FindBin(ts+fDeltaTVD);
1272 0 : int ntr = ntrCumul[bmaxT]-ntrCumul[bminT];
1273 0 : int binv = (ts+1-fTMin)/fDeltaTVD;
1274 0 : nChecked++;
1275 0 : if (ntr>fMinTracksPerVDBin) {
1276 0 : binStat[binv] = kTRUE; // flag complete bin
1277 0 : if (prevEmpty>longestEmpty) longestEmpty = prevEmpty;
1278 : prevEmpty = 0; // reset empty segment counter
1279 0 : continue;
1280 : }
1281 : // do we expect enough tracks in this time period?
1282 0 : int nexpTracks = (nevCumul[bmaxE]-nevCumul[bminE])*fEstTracksPerEvent;
1283 : // declare stat insufficient only if we expect enough tracks there
1284 0 : if (nexpTracks>fMinTracksPerVDBin && nexpTracks>nexpTracksAv/2 && ntr<fMinTracksPerVDBin/2) {
1285 0 : nHoles++;
1286 0 : prevEmpty += fDeltaTVD;
1287 0 : }
1288 0 : else binStat[binv] = kTRUE; // flag complete bin, since this bin is depleted on event level
1289 0 : }
1290 0 : if (prevEmpty>longestEmpty) longestEmpty = prevEmpty;
1291 0 : AliInfoF("%d bins out %d checked got enough stat., longest empty segment: %ds",
1292 : nChecked-nHoles,nChecked,longestEmpty);
1293 0 : if (nChecked && nHoles<maxHolesFrac*nChecked && longestEmpty<fSigmaTVD*0.75) return kTRUE;
1294 0 : return kFALSE;
1295 0 : }
1296 :
1297 : //________________________________________________
1298 : void AliTPCDcalibRes::FillDriftResidualsTrees()
1299 : {
1300 : // fill local trees for vdrift calibration
1301 0 : fDTV.t = fTimeStamp;
1302 0 : for (int icl=fNCl;icl--;) {
1303 0 : if (fArrR[icl]<kInvalidR) continue; // rejected outlier
1304 0 : Bool_t isCside = ((fArrSectID[icl]/kNSect)&0x1);
1305 0 : fDTV.side = isCside ? -1:1;
1306 0 : fDTV.dz = fArrDZ[icl];
1307 0 : fDTV.drift = kZLim[isCside] - fDTV.side*fArrZCl[icl];
1308 0 : fDTV.ylab = fArrR[icl]*TMath::Sin(fArrPhi[icl]);
1309 0 : fDTV.r = fArrR[icl];
1310 : //
1311 0 : fTmpTree[0]->Fill();
1312 : }
1313 0 : if (fTracksRate) fTracksRate->Fill(fTimeStamp); // register track time
1314 : //
1315 0 : }
1316 :
1317 : //________________________________________________
1318 : void AliTPCDcalibRes::FillLocalResidualsTrees()
1319 : {
1320 : // fill local trees with binned data
1321 0 : float voxVars[kVoxHDim]={0}; // voxel variables (unbinned)
1322 0 : for (int icl=fNCl;icl--;) {
1323 0 : if (fArrX[icl]<kInvalidR) continue; // rejected outlier
1324 0 : int sectID = fArrSectID[icl]; // 0-35 numbering
1325 : //
1326 : // calculate voxel variables and bins
1327 : //
1328 0 : if (!FindVoxelBin(sectID, fArrX[icl], fArrYCl[icl], fArrZCl[icl], fDTS.bvox, voxVars)) continue;
1329 0 : fDTS.dy = fArrDY[icl];
1330 0 : fDTS.dz = fArrDZ[icl];
1331 0 : fDTS.tgSlp = fArrTgSlp[icl];
1332 : //
1333 0 : fTmpTree[sectID]->Fill();
1334 : //
1335 : // fill statistics on distribution within the voxel, last dimension, kVoxV is for Nentries
1336 0 : ULong64_t binToFill = GetBin2Fill(fDTS.bvox,kVoxV); // bin of sector stat histo
1337 0 : float &binEntries = fArrNDStat[sectID]->At(binToFill); // entries in the voxel
1338 0 : float oldEntries = binEntries++;
1339 0 : float norm = 1.f/binEntries;
1340 0 : for (int iv=kVoxDim;iv--;) {
1341 0 : float &mean = fArrNDStat[sectID]->At(binToFill+iv-kVoxV);
1342 0 : mean = ( mean*oldEntries + voxVars[iv]) * norm; // account new bin entry in averages calculation
1343 : }
1344 : //
1345 0 : } // loop over clusters
1346 : //
1347 0 : if (fTracksRate) fTracksRate->Fill(fTimeStamp); // register track time
1348 : //
1349 0 : }
1350 :
1351 : //________________________________________________
1352 : void AliTPCDcalibRes::FillCorrectedResiduals()
1353 : {
1354 : // fill local trees result of closure test: corrected distortions
1355 :
1356 0 : float voxVars[kVoxHDim]={0}; // voxel variables (unbinned)
1357 0 : fDTC.t = fTimeStamp;
1358 0 : fDTC.q2pt = fQ2Pt;
1359 0 : fDTC.tgLam = fTgLam;
1360 : //
1361 0 : for (int icl=fNCl;icl--;) {
1362 0 : if (fArrX[icl]<kInvalidR) continue; // rejected outlier
1363 0 : int sectID = fArrSectID[icl]; // 0-35 numbering
1364 : //
1365 : // extract correction
1366 : // calculate voxel variables and bins
1367 0 : if (!FindVoxelBin(sectID,fArrX[icl], fArrYCl[icl], fArrZCl[icl], fDTC.bvox, voxVars)) continue;
1368 0 : int row159 = GetRowID(fArrX[icl]);
1369 0 : if (row159<0) continue;
1370 0 : float corr[3];
1371 :
1372 0 : fChebCorr->Eval(sectID, row159, fArrYCl[icl]/fArrX[icl], fArrZCl[icl]/fArrX[icl], corr);
1373 : //
1374 0 : fDTC.dyR = fArrDY[icl];
1375 0 : fDTC.dzR = fArrDZ[icl];
1376 :
1377 0 : fDTC.dyC = fArrDY[icl] - (corr[kResY]-corr[kResX]*fArrTgSlp[icl]);
1378 0 : fDTC.dzC = fArrDZ[icl] - (corr[kResZ]-corr[kResX]*fTgLam); // we evaluate at pad-row
1379 :
1380 0 : fDTC.tgSlp = fArrTgSlp[icl];
1381 0 : fDTC.x = fArrX[icl];
1382 0 : fDTC.y = fArrYCl[icl];
1383 0 : fDTC.z = fArrZCl[icl];
1384 : //
1385 0 : fTmpTree[sectID]->Fill();
1386 : //
1387 0 : } // loop over clusters
1388 0 : }
1389 :
1390 : //________________________________________________
1391 : void AliTPCDcalibRes::CreateLocalResidualsTrees(int mode)
1392 : {
1393 : // temporary trees for local delta's storage
1394 : //
1395 0 : static dts_t *dtsP = &fDTS;
1396 0 : static dtc_t *dtcP = &fDTC;
1397 0 : static dtv_t *dtvP = &fDTV;
1398 0 : TString namef;
1399 0 : if (mode==kVDriftCalibMode) {
1400 0 : namef = Form("%s.root",kDriftResFileName);
1401 0 : fTmpFile[0] = TFile::Open(namef.Data(),"recreate");
1402 0 : fTmpTree[0] = new TTree("resdrift","");
1403 0 : fTmpTree[0]->Branch("dtv", &dtvP);
1404 : }
1405 0 : else if (mode==kDistExtractMode||mode==kDistClosureTestMode) {
1406 0 : for (int is=0;is<kNSect2;is++) {
1407 0 : if (mode==kDistExtractMode) namef = Form("%s%d.root",kLocalResFileName,is);
1408 0 : else /*if (mode==kDistClosureTestMode)*/ namef = Form("%s%d.root",kClosureTestFileName,is);
1409 0 : fTmpFile[is] = TFile::Open(namef.Data(),"recreate");
1410 0 : fTmpTree[is] = new TTree(Form("ts%d",is),"");
1411 : //
1412 0 : if (mode==kDistExtractMode) {
1413 0 : fTmpTree[is]->Branch("dts", &dtsP);
1414 : //fTmpTree[is]->SetAutoFlush(150000);
1415 : //
1416 0 : fStatHist[is] = CreateVoxelStatHisto(is);
1417 0 : fArrNDStat[is] = (TNDArrayT<float>*)&fStatHist[is]->GetArray();
1418 0 : }
1419 0 : else if (mode==kDistClosureTestMode) {
1420 0 : fTmpTree[is]->Branch("dtc", &dtcP);
1421 : }
1422 : }
1423 0 : }
1424 0 : else AliFatalF("Unknown mode %d",mode);
1425 : //
1426 0 : }
1427 :
1428 : //________________________________________________
1429 : void AliTPCDcalibRes::CloseLocalResidualsTrees(int /*mode*/)
1430 : {
1431 : // close trees for local delta's storage
1432 : //
1433 0 : for (int is=0;is<kNSect2;is++) {
1434 0 : if (!fTmpFile[is]) continue;
1435 0 : fTmpFile[is]->cd();
1436 0 : fTmpTree[is]->Write("", TObject::kOverwrite);
1437 0 : delete fTmpTree[is];
1438 0 : fTmpTree[is] = 0;
1439 0 : fTmpFile[is]->Close();
1440 0 : delete fTmpFile[is];
1441 0 : fTmpFile[is] = 0;
1442 0 : }
1443 : //
1444 0 : }
1445 :
1446 : //__________________________________________________________________________________
1447 : Bool_t AliTPCDcalibRes::CompareToHelix(float *resHelixY, float *resHelixZ)
1448 : {
1449 : // compare track to helix, refit q/pt and tgLambda and build array of tg(slope) at pad-rows
1450 : const double kEps = 1e-12;
1451 0 : float xlab[kNPadRows],ylab[kNPadRows],spath[kNPadRows]; // lab X,Y rotated to for sectort of 1st cluster
1452 : // fill lab coordinates
1453 0 : float crv = TMath::Abs(fQ2Pt*fBz*0.299792458e-3f), cs,sn;
1454 0 : int sectPrev=-1,sect0 = fArrSectID[0]%kNSect; // align to the sector of 1st point
1455 0 : float phiSect = (sect0+0.5)*20*TMath::DegToRad();
1456 0 : double sna = TMath::Sin(phiSect), csa = TMath::Cos(phiSect);
1457 : //
1458 0 : spath[0] = 0.f;
1459 0 : for (int ip=0;ip<fNCl;ip++) {
1460 0 : cs = TMath::Cos(fArrPhi[ip]-phiSect);
1461 0 : sn = TMath::Sin(fArrPhi[ip]-phiSect);
1462 0 : xlab[ip] = fArrR[ip]*cs - fArrDY[ip]*sn;
1463 0 : ylab[ip] = fArrDY[ip]*cs + fArrR[ip]*sn;
1464 0 : if (ip) {
1465 0 : float dx = xlab[ip]-xlab[ip-1];
1466 0 : float dy = ylab[ip]-ylab[ip-1];
1467 0 : float ds2 = dx*dx+dy*dy;
1468 0 : float ds = TMath::Sqrt(ds2); // circular path
1469 0 : if (ds*crv>0.05) {
1470 : // account for the arc-chord difference as 1st 2 terms of asin expansion
1471 0 : ds *= (1.f+ds2*crv*crv/24.f);
1472 0 : }
1473 0 : spath[ip] = spath[ip-1]+ds;
1474 0 : }
1475 : }
1476 0 : double xcSec=0,ycSec=0,xc=0,yc=0,r=0;
1477 0 : FitCircle(fNCl,xlab,ylab,xcSec,ycSec,r,resHelixY);
1478 : // determine qurvature
1479 0 : float phi0 = TMath::ATan2(ylab[0],xlab[0]);
1480 0 : if (phi0<0) phi0 += TMath::Pi()*2;
1481 0 : float phi1 = TMath::ATan2(ylab[fNCl-1],xlab[fNCl-1]);
1482 0 : if (phi1<0) phi1 += TMath::Pi()*2;
1483 0 : float dphi = phi1-phi0;
1484 : int curvSign = 1;
1485 0 : if (dphi>0) {
1486 0 : if (dphi<TMath::Pi()) curvSign = -1; // clockwise, no 2pi-0 crossing
1487 : }
1488 0 : else if (dphi<-TMath::Pi()) curvSign = -1; // clockwise, 2pi-0 crossing
1489 : //
1490 0 : fQ2Pt = curvSign/(r*fBz*0.299792458e-3f);
1491 : //
1492 : // calculate circle coordinates in the lab frame
1493 0 : xc = xcSec*csa - ycSec*sna;
1494 0 : yc = ycSec*csa + xcSec*sna;
1495 : //
1496 0 : float pol1z[2],pol1zE[4] ;
1497 0 : Bool_t resfZ = FitPoly1(spath, fArrZTr, 0, fNCl, pol1z, pol1zE);
1498 : //
1499 0 : fTgLam = pol1z[1]; // new tg. lambda
1500 : // extract deviations wrt helical fit and fill track slopes in sector frame
1501 : float hmnY=1e9,hmxY=-1e9,hmnZ=1e9,hmxZ=-1e9;
1502 :
1503 0 : for (int ip=0;ip<fNCl;ip++) {
1504 0 : float val = fArrZTr[ip] - (pol1z[0]+spath[ip]*pol1z[1]);
1505 0 : resHelixZ[ip] = val;
1506 0 : if (val<hmnZ) hmnZ = val;
1507 0 : if (val>hmxZ) hmxZ = val;
1508 : //
1509 0 : val = resHelixY[ip];
1510 0 : if (val<hmnY) hmnY = val;
1511 0 : if (val>hmxY) hmxY = val;
1512 : //
1513 0 : int sect = fArrSectID[ip]%kNSect;
1514 0 : if (sect!=sect0) {
1515 : sect0 = sect;
1516 0 : phiSect = (0.5f + sect)*kSecDPhi;
1517 0 : sna = TMath::Sin(phiSect);
1518 0 : csa = TMath::Cos(phiSect);
1519 0 : xcSec = xc*csa + yc*sna; // recalculate circle center in the sector frame
1520 0 : }
1521 : // find intersection of the circle with the padrow
1522 : // 1) equation of circle in lab: x=xc+r*cos(tau), y=yc+r*sin(tau)
1523 : // 2) equation of circle in sector frame:
1524 : // x=xc'+R*cos(tau-alpSect), y=yc'+R*sin(tau-alpSect)
1525 : // with xc'=xc*cos(-alp)-yc*sin(-alp); yc'=yc*cos(-alp)+xc*sin(-alp)
1526 : // The circle and padrow at X cross at cos(tau) = (X-xc*csa+yc*sna)/R
1527 : // Hence the derivative of y vs x in sector frame:
1528 0 : cs = TMath::Cos(fArrPhi[ip]-phiSect);
1529 0 : double xRow = fArrR[ip]*cs;
1530 0 : double cstalp = (xRow - xcSec)/r;
1531 0 : if (TMath::Abs(cstalp)>1.-kEps) { // track cannot reach this padrow
1532 0 : cstalp = TMath::Sign(1.-kEps,cstalp);
1533 0 : }
1534 : // and the slope in sector frame is +-1/tg(acos(cstalp)) = +-cstalp/sqrt(1-cstalp^2)
1535 : // The sign is defined by the fact that in B+ the slope of q- should increase with X.
1536 : // Since the derivative of cstalp/sqrt(1-cstalp^2) on X is positive, just look on qB
1537 0 : fArrTgSlp[ip] = cstalp/TMath::Sqrt((1.-cstalp)*(1.+cstalp));
1538 0 : if (fQ2Pt*fBz>0) fArrTgSlp[ip] = -fArrTgSlp[ip];
1539 : }
1540 : //
1541 : // if (TMath::Abs(hmxY-hmnY)>fMaxDevYHelix || TMath::Abs(hmxZ-hmnZ)>fMaxDevZHelix)
1542 : // printf("MinMax%d: %e %e %e %e\n",evID,hmnY,hmxY,hmnZ,hmxZ);
1543 0 : return TMath::Abs(hmxY-hmnY)<fMaxDevYHelix && TMath::Abs(hmxZ-hmnZ)<fMaxDevZHelix;
1544 0 : }
1545 :
1546 : //________________________________________________
1547 : void AliTPCDcalibRes::ClosureTest()
1548 : {
1549 : // correct distortions
1550 0 : TStopwatch sw;
1551 0 : sw.Start();
1552 0 : if (!fChebCorr) {
1553 0 : AliError("Chebyshev correction object was not created, cannot run closure test");
1554 0 : return;
1555 : }
1556 0 : CollectData(kDistClosureTestMode);
1557 0 : sw.Stop();
1558 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
1559 : //
1560 0 : }
1561 :
1562 : //________________________________________________
1563 : void AliTPCDcalibRes::ProcessResiduals()
1564 : {
1565 : // project local trees, extract distortions
1566 0 : if (!fInitDone) Init(); //{AliError("Init not done"); return;}
1567 0 : LoadStatHistos();
1568 0 : AliSysInfo::AddStamp("ProcResid",0,0,0,0);
1569 : //
1570 0 : for (int is=0;is<kNSect2;is++) ProcessSectorResiduals(is);
1571 : //
1572 0 : AliSysInfo::AddStamp("ProcResid",1,0,0,0);
1573 : //
1574 0 : }
1575 :
1576 : //________________________________________________
1577 : void AliTPCDcalibRes::ProcessDispersions()
1578 : {
1579 : // extract distortions of corrected Y residuals ||| DEPRECATED
1580 0 : if (!fInitDone) Init(); //{AliError("Init not done"); return;}
1581 : //
1582 0 : LoadStatHistos();
1583 0 : for (int is=0;is<kNSect2;is++) {
1584 0 : ProcessSectorDispersions(is);
1585 0 : if (fDeleteSectorTrees) {
1586 0 : TString sectFileName = Form("%s%d.root",kLocalResFileName,is);
1587 0 : AliInfoF("Deleting %s",sectFileName.Data());
1588 0 : unlink(sectFileName.Data());
1589 0 : }
1590 : }
1591 : //
1592 0 : }
1593 :
1594 : //________________________________________________
1595 : void AliTPCDcalibRes::ProcessSectorDispersions(int is)
1596 : {
1597 : // extract dispersion of corrected residuals ||| DEPRECATED
1598 : const float kEps = 1e-6;
1599 :
1600 0 : if (!fInitDone) {AliError("Init not done"); return;}
1601 0 : TStopwatch sw; sw.Start();
1602 0 : AliSysInfo::AddStamp("ProcessSectorDispersions",is,0,0,0);
1603 :
1604 0 : TString sectFileName = Form("%s%d.root",kLocalResFileName,is);
1605 0 : TFile* sectFile = TFile::Open(sectFileName.Data());
1606 0 : if (!sectFile) AliFatalF("file %s not found",sectFileName.Data());
1607 0 : TString treeName = Form("ts%d",is);
1608 0 : TTree *sectTree = (TTree*) sectFile->Get(treeName.Data());
1609 0 : if (!sectTree) AliFatalF("tree %s is not found in file %s",treeName.Data(),sectFileName.Data());
1610 : //
1611 0 : dts_t *dtsP = &fDTS;
1612 0 : sectTree->SetBranchAddress("dts",&dtsP);
1613 0 : int npoints = sectTree->GetEntries();
1614 0 : if (!npoints) {
1615 0 : AliWarningF("No entries for sector %d",is);
1616 0 : delete sectTree;
1617 0 : sectFile->Close(); // to reconsider: reuse the file
1618 0 : delete sectFile;
1619 0 : return;
1620 : }
1621 0 : Short_t *resYArr = new Short_t[npoints];
1622 0 : Short_t *tgslArr = new Short_t[npoints];
1623 0 : UShort_t *binArr = new UShort_t[npoints];
1624 0 : Int_t* index = new Int_t[npoints];
1625 0 : TArrayF dya(1000),tga(1000);//, dza(1000),
1626 0 : float *dy = dya.GetArray(), *tg = tga.GetArray();//, *dz = dza.GetArray();
1627 : int nacc = 0;
1628 0 : bres_t* sectData = fSectGVoxRes[is];
1629 0 : for (int ie=0;ie<npoints;ie++) {
1630 0 : sectTree->GetEntry(ie);
1631 0 : if (TMath::Abs(fDTS.tgSlp)>=kMaxTgSlp) continue;
1632 0 : resYArr[nacc] = Short_t(fDTS.dy*0x7fff/kMaxResid);
1633 0 : tgslArr[nacc] = Short_t(fDTS.tgSlp*0x7fff/kMaxTgSlp);
1634 0 : binArr[nacc] = GetVoxGBin(fDTS.bvox);
1635 0 : nacc++;
1636 0 : }
1637 0 : TMath::Sort(nacc, binArr, index, kFALSE); // sort in voxel increasing order
1638 : UShort_t curBin = 0xffff;
1639 0 : UChar_t bvox[kVoxDim];
1640 : int nproc = 0, npBin = 0;
1641 0 : while (nproc<nacc) {
1642 0 : int ip = index[nproc++];
1643 0 : if (curBin!=binArr[ip]) {
1644 0 : if (npBin) {
1645 0 : bres_t& resVox = sectData[curBin];
1646 0 : GBin2Vox(curBin,resVox.bvox); // parse voxel
1647 0 : ProcessVoxelDispersions(npBin,tg,dy,resVox);
1648 0 : }
1649 0 : curBin = binArr[ip];
1650 : npBin = 0;
1651 0 : }
1652 0 : if (npBin==dya.GetSize()) {
1653 0 : dya.Set(100+npBin); dy = dya.GetArray();
1654 0 : tga.Set(100+npBin); tg = tga.GetArray();
1655 0 : }
1656 0 : dy[npBin] = resYArr[ip]*kMaxResid/0x7fff;
1657 0 : tg[npBin] = tgslArr[ip]*kMaxTgSlp/0x7fff;
1658 0 : npBin++;
1659 : }
1660 0 : if (npBin) {
1661 0 : bres_t& resVox = sectData[curBin];
1662 0 : GBin2Vox(curBin,resVox.bvox); // parse voxel
1663 0 : ProcessVoxelDispersions(npBin,tg,dy,resVox);
1664 0 : }
1665 : //
1666 0 : delete[] binArr;
1667 0 : delete[] resYArr;
1668 0 : delete[] tgslArr;
1669 0 : delete[] index;
1670 : //
1671 0 : delete sectTree;
1672 0 : sectFile->Close(); // to reconsider: reuse the file
1673 0 : delete sectFile;
1674 : //
1675 : // now smooth the dispersion
1676 0 : for (bvox[kVoxX]=0;bvox[kVoxX]<fNXBins;bvox[kVoxX]++) {
1677 0 : if (GetXBinIgnored(is,bvox[kVoxX])) continue;
1678 0 : for (bvox[kVoxZ]=0;bvox[kVoxZ]<fNZ2XBins;bvox[kVoxZ]++) {
1679 0 : for (bvox[kVoxF]=0;bvox[kVoxF]<fNY2XBins;bvox[kVoxF]++) {
1680 0 : int binGlo = GetVoxGBin(bvox);
1681 0 : bres_t *voxRes = §Data[binGlo];
1682 0 : Bool_t res = GetSmoothEstimate(is,voxRes->stat[kVoxX],voxRes->stat[kVoxF],voxRes->stat[kVoxZ],
1683 0 : BIT(kResD), voxRes->DS);
1684 : }
1685 : }
1686 : }
1687 : //
1688 0 : sw.Stop();
1689 0 : AliInfoF("Sector%2d | timing: real: %.3f cpu: %.3f",is, sw.RealTime(), sw.CpuTime());
1690 0 : AliSysInfo::AddStamp("ProcessSectorDispersions",1,0,0,0);
1691 : //
1692 0 : }
1693 :
1694 : //_________________________________________________
1695 : void AliTPCDcalibRes::ProcessSectorResiduals(int is)
1696 : {
1697 : // process residuals for single sector staring from local binned per-sector trees
1698 : //
1699 : const int kMaxPnt = 30000000; // max points per sector to accept
1700 0 : TStopwatch sw; sw.Start();
1701 0 : AliSysInfo::AddStamp("ProcSectRes",is,0,0,0);
1702 : //
1703 0 : fNSmoothingFailedBins[is] = 0;
1704 : //
1705 0 : TString sectFileName = Form("%s%d.root",kLocalResFileName,is);
1706 0 : TFile* sectFile = TFile::Open(sectFileName.Data());
1707 0 : if (!sectFile) AliFatalF("file %s not found",sectFileName.Data());
1708 0 : TString treeName = Form("ts%d",is);
1709 0 : TTree *sectTree = (TTree*) sectFile->Get(treeName.Data());
1710 0 : if (!sectTree) AliFatalF("tree %s is not found in file %s",treeName.Data(),sectFileName.Data());
1711 : //
1712 0 : if (fSectGVoxRes[is]) delete[] fSectGVoxRes[is];
1713 0 : fSectGVoxRes[is] = new bres_t[fNGVoxPerSector]; // here we keep main result
1714 : bres_t* sectData = fSectGVoxRes[is];
1715 : // by default set the COG estimates to bin center
1716 0 : for (int ix=0;ix<fNXBins;ix++) {
1717 0 : for (int ip=0;ip<fNY2XBins;ip++) {
1718 0 : for (int iz=0;iz<fNZ2XBins;iz++) { // extract line in z
1719 0 : int binGlo = GetVoxGBin(ix,ip,iz);
1720 0 : bres_t &resVox = sectData[binGlo];
1721 0 : resVox.bvox[kVoxX] = ix;
1722 0 : resVox.bvox[kVoxF] = ip;
1723 0 : resVox.bvox[kVoxZ] = iz;
1724 0 : resVox.bsec = is;
1725 0 : GetVoxelCoordinates(resVox.bsec,resVox.bvox[kVoxX],resVox.bvox[kVoxF],resVox.bvox[kVoxZ],
1726 0 : resVox.stat[kVoxX],resVox.stat[kVoxF],resVox.stat[kVoxZ]);
1727 : }
1728 : }
1729 : }
1730 : //
1731 0 : dts_t *dtsP = &fDTS;
1732 0 : sectTree->SetBranchAddress("dts",&dtsP);
1733 0 : int npoints = sectTree->GetEntries();
1734 0 : if (!npoints) {
1735 0 : AliWarningF("No entries for sector %d, masking all rows",is);
1736 0 : for (int ix=fNXBins;ix--;) SetXBinIgnored(is,ix);
1737 0 : delete sectTree;
1738 0 : sectFile->Close(); // to reconsider: reuse the file
1739 0 : delete sectFile;
1740 0 : return;
1741 : }
1742 0 : if (npoints>kMaxPnt) npoints = kMaxPnt;
1743 0 : sw.Stop();
1744 0 : AliInfoF("Sector%2d. Extracted %d points of unbinned data. Timing: real: %.3f cpu: %.3f",
1745 : is, npoints, sw.RealTime(), sw.CpuTime());
1746 0 : sw.Start(kFALSE);
1747 : //
1748 0 : Short_t *resYArr = new Short_t[npoints];
1749 0 : Short_t *resZArr = new Short_t[npoints];
1750 0 : Short_t *tgslArr = new Short_t[npoints];
1751 0 : UShort_t *binArr = new UShort_t[npoints];
1752 0 : Int_t* index = new Int_t[npoints];
1753 0 : TArrayF dya(1000),dza(1000),tga(1000);
1754 0 : float *dy = dya.GetArray(), *dz = dza.GetArray(), *tg = tga.GetArray();
1755 : int nacc = 0;
1756 0 : for (int ie=0;ie<npoints;ie++) {
1757 0 : sectTree->GetEntry(ie);
1758 0 : if (TMath::Abs(fDTS.tgSlp)>=kMaxTgSlp) continue;
1759 0 : resYArr[nacc] = Short_t(fDTS.dy*0x7fff/kMaxResid);
1760 0 : resZArr[nacc] = Short_t(fDTS.dz*0x7fff/kMaxResid);
1761 0 : tgslArr[nacc] = Short_t(fDTS.tgSlp*0x7fff/kMaxTgSlp);
1762 0 : binArr[nacc] = GetVoxGBin(fDTS.bvox);
1763 0 : nacc++;
1764 0 : }
1765 : //
1766 0 : delete sectTree;
1767 0 : sectFile->Close(); // to reconsider: reuse the file
1768 0 : delete sectFile;
1769 : //
1770 0 : TMath::Sort(nacc, binArr, index, kFALSE); // sort in voxel increasing order
1771 : UShort_t curBin = 0xffff;
1772 0 : UChar_t bvox[kVoxDim];
1773 : int nproc = 0, npBin = 0;
1774 0 : while (nproc<nacc) {
1775 0 : int ip = index[nproc++];
1776 0 : if (curBin!=binArr[ip]) {
1777 0 : if (npBin) {
1778 0 : bres_t& resVox = sectData[curBin];
1779 0 : GBin2Vox(curBin,resVox.bvox); // parse voxel
1780 0 : ProcessVoxelResiduals(npBin,tg,dy,dz,resVox);
1781 0 : }
1782 0 : curBin = binArr[ip];
1783 : npBin = 0;
1784 0 : }
1785 0 : if (npBin==dya.GetSize()) {
1786 0 : dya.Set(100+npBin); dy = dya.GetArray();
1787 0 : dza.Set(100+npBin); dz = dza.GetArray();
1788 0 : tga.Set(100+npBin); tg = tga.GetArray();
1789 0 : }
1790 0 : dy[npBin] = resYArr[ip]*kMaxResid/0x7fff;
1791 0 : dz[npBin] = resZArr[ip]*kMaxResid/0x7fff;
1792 0 : tg[npBin] = tgslArr[ip]*kMaxTgSlp/0x7fff;
1793 0 : npBin++;
1794 : }
1795 0 : if (npBin) {
1796 0 : bres_t &resVox = sectData[curBin];
1797 0 : GBin2Vox(curBin,resVox.bvox); // parse voxel
1798 0 : ProcessVoxelResiduals(npBin,tg,dy,dz,resVox);
1799 0 : }
1800 : //
1801 0 : sw.Stop();
1802 0 : AliInfoF("Sector%2d. Extracted residuals. Timing: real: %.3f cpu: %.3f",
1803 : is, sw.RealTime(), sw.CpuTime());
1804 0 : sw.Start(kFALSE);
1805 :
1806 0 : int nrowOK = ValidateVoxels(is);
1807 0 : if (!nrowOK) AliWarningF("Sector%2d: all X-bins disabled, abandon smoothing",is);
1808 0 : else Smooth0(is);
1809 : //
1810 0 : sw.Stop();
1811 0 : AliInfoF("Sector%2d. Smoothed residuals. Timing: real: %.3f cpu: %.3f",
1812 : is, sw.RealTime(), sw.CpuTime());
1813 0 : sw.Start(kFALSE);
1814 :
1815 : // now process dispersions
1816 : curBin = 0xffff;
1817 :
1818 : nproc = 0;
1819 : npBin = 0;
1820 0 : while (nproc<nacc) {
1821 0 : int ip = index[nproc++];
1822 0 : if (curBin!=binArr[ip]) {
1823 0 : if (npBin) {
1824 0 : bres_t& resVox = sectData[curBin];
1825 0 : GBin2Vox(curBin,resVox.bvox); // parse voxel
1826 0 : if (!GetXBinIgnored(is,resVox.bvox[kVoxX])) ProcessVoxelDispersions(npBin,tg,dy,resVox);
1827 0 : }
1828 0 : curBin = binArr[ip];
1829 : npBin = 0;
1830 0 : }
1831 0 : if (npBin==dya.GetSize()) {
1832 0 : dya.Set(100+npBin); dy = dya.GetArray();
1833 0 : tga.Set(100+npBin); tg = tga.GetArray();
1834 0 : }
1835 0 : dy[npBin] = resYArr[ip]*kMaxResid/0x7fff;
1836 0 : tg[npBin] = tgslArr[ip]*kMaxTgSlp/0x7fff;
1837 0 : npBin++;
1838 : }
1839 0 : if (npBin) {
1840 0 : bres_t& resVox = sectData[curBin];
1841 0 : GBin2Vox(curBin,resVox.bvox); // parse voxel
1842 0 : if (!GetXBinIgnored(is,resVox.bvox[kVoxX])) ProcessVoxelDispersions(npBin,tg,dy,resVox);
1843 0 : }
1844 : //
1845 : // now smooth the dispersion
1846 0 : for (bvox[kVoxX]=0;bvox[kVoxX]<fNXBins;bvox[kVoxX]++) {
1847 0 : if (GetXBinIgnored(is,bvox[kVoxX])) continue;
1848 0 : for (bvox[kVoxZ]=0;bvox[kVoxZ]<fNZ2XBins;bvox[kVoxZ]++) {
1849 0 : for (bvox[kVoxF]=0;bvox[kVoxF]<fNY2XBins;bvox[kVoxF]++) {
1850 0 : int binGlo = GetVoxGBin(bvox);
1851 0 : bres_t *voxRes = §Data[binGlo];
1852 0 : Bool_t res = GetSmoothEstimate(is,voxRes->stat[kVoxX],voxRes->stat[kVoxF],voxRes->stat[kVoxZ],
1853 0 : BIT(kResD), voxRes->DS);
1854 : }
1855 : }
1856 : }
1857 :
1858 0 : delete[] binArr;
1859 0 : delete[] resYArr;
1860 0 : delete[] resZArr;
1861 0 : delete[] tgslArr;
1862 0 : delete[] index;
1863 : //
1864 0 : sw.Stop();
1865 0 : AliInfoF("Sector%2d. Processed dispersion. Timing: real: %.3f cpu: %.3f",is, sw.RealTime(), sw.CpuTime());
1866 0 : AliSysInfo::AddStamp("ProcSectRes",is,1,0,0);
1867 : //
1868 0 : }
1869 : //________________________________________________
1870 : void AliTPCDcalibRes::ReProcessResiduals()
1871 : {
1872 : // reprocess residuals using raw voxel info filled from existing resVoxTree
1873 : // The raw data is already loaded from the tree
1874 0 : AliSysInfo::AddStamp("ReProcResid",0,0,0,0);
1875 0 : for (int is=0;is<kNSect2;is++) ReProcessSectorResiduals(is);
1876 0 : AliSysInfo::AddStamp("ReProcResid",1,0,0,0);
1877 0 : }
1878 :
1879 : //_________________________________________________
1880 : void AliTPCDcalibRes::ReProcessSectorResiduals(int is)
1881 : {
1882 : // Reprocess residuals for single sector filled from existing resVoxTree
1883 : // The raw data is already loaded from the tree
1884 : //
1885 0 : TStopwatch sw; sw.Start();
1886 0 : AliSysInfo::AddStamp("RProcSectRes",is,0,0,0);
1887 : //
1888 0 : fNSmoothingFailedBins[is] = 0;
1889 0 : bres_t* sectData = fSectGVoxRes[is];
1890 0 : if (!sectData) AliFatalF("No SectGVoxRes data for sector %d",is);
1891 : //
1892 0 : int nrowOK = ValidateVoxels(is);
1893 0 : if (!nrowOK) AliWarningF("Sector%2d: all X-bins disabled, abandon smoothing",is);
1894 0 : else Smooth0(is);
1895 : //
1896 0 : UChar_t bvox[kVoxDim];
1897 : // now smooth the dispersion
1898 0 : for (bvox[kVoxX]=0;bvox[kVoxX]<fNXBins;bvox[kVoxX]++) {
1899 0 : if (GetXBinIgnored(is,bvox[kVoxX])) continue;
1900 0 : for (bvox[kVoxZ]=0;bvox[kVoxZ]<fNZ2XBins;bvox[kVoxZ]++) {
1901 0 : for (bvox[kVoxF]=0;bvox[kVoxF]<fNY2XBins;bvox[kVoxF]++) {
1902 0 : int binGlo = GetVoxGBin(bvox);
1903 0 : bres_t *voxRes = §Data[binGlo];
1904 0 : Bool_t res = GetSmoothEstimate(is,voxRes->stat[kVoxX],voxRes->stat[kVoxF],voxRes->stat[kVoxZ],
1905 0 : BIT(kResD), voxRes->DS);
1906 : }
1907 : }
1908 : }
1909 0 : sw.Stop();
1910 0 : AliInfoF("Sector%2d. Processed dispersion. Timing: real: %.3f cpu: %.3f",is, sw.RealTime(), sw.CpuTime());
1911 0 : AliSysInfo::AddStamp("ReProcSectRes",is,1,0,0);
1912 : //
1913 0 : }
1914 :
1915 : //_________________________________________________________
1916 : Float_t AliTPCDcalibRes::FitPoly1Robust(int np, float* x, float* y, float* res, float* err, float ltmCut)
1917 : {
1918 : // robust pol1 fit, modifies input arrays order
1919 0 : res[0] = res[1] = 0.f;
1920 0 : if (np<2) return -1;
1921 0 : TVectorF yres(7);
1922 0 : int *indY = TStatToolkit::LTMUnbinned(np,y,yres,ltmCut);
1923 0 : if (!indY) return -1;
1924 : // rearrange used events in increasing order
1925 0 : TStatToolkit::Reorder(np,y,indY);
1926 0 : TStatToolkit::Reorder(np,x,indY);
1927 : //
1928 : // 1st fit to get crude slope
1929 0 : int npuse = TMath::Nint(yres[0]);
1930 0 : int offs = TMath::Nint(yres[5]);
1931 : // use only entries selected by LTM for the fit
1932 0 : float a,b;
1933 0 : AliTPCDcalibRes::medFit(npuse, x+offs, y+offs, a, b, err);
1934 : //
1935 : // don't abuse stack
1936 0 : float *ycmHeap=0,ycmStack[np<kMaxOnStack ? np:1],*ycm=np<kMaxOnStack ? &ycmStack[0] : (ycmHeap=new float[np]);
1937 0 : int *indcmHeap=0,indcmStack[np<kMaxOnStack ? np:1],*indcm=np<kMaxOnStack ? &indcmStack[0] : (indcmHeap=new int[np]);
1938 : //
1939 0 : for (int i=np;i--;) ycm[i] = y[i]-(a+b*x[i]);
1940 0 : TMath::Sort(np,ycm,indcm,kFALSE);
1941 0 : TStatToolkit::Reorder(np,ycm,indcm);
1942 0 : TStatToolkit::Reorder(np,y,indcm); // we must keep the same order
1943 0 : TStatToolkit::Reorder(np,x,indcm);
1944 : //
1945 : // robust estimate of sigma after crude slope correction
1946 0 : float sigMAD = AliTPCDcalibRes::MAD2Sigma(npuse,ycm+offs);
1947 : // find LTM estimate matching to sigMAD, keaping at least given fraction
1948 0 : indY = AliTPCDcalibRes::LTMUnbinnedSig(np, ycm, yres, sigMAD,0.5,kTRUE);
1949 0 : delete[] ycmHeap;
1950 0 : delete[] indcmHeap;
1951 : //
1952 0 : if (!indY) return -1;
1953 : // final fit
1954 0 : npuse = TMath::Nint(yres[0]);
1955 0 : offs = TMath::Nint(yres[5]);
1956 0 : AliTPCDcalibRes::medFit(npuse, x+offs, y+offs, a,b, err);
1957 0 : res[0] = a;
1958 0 : res[1] = b;
1959 0 : return sigMAD;
1960 0 : }
1961 :
1962 : //_________________________________________________
1963 : void AliTPCDcalibRes::ProcessVoxelResiduals(int np, float* tg, float *dy, float *dz, bres_t& voxRes)
1964 : {
1965 : // extract X,Y,Z distortions of the voxel
1966 0 : if (np<fMinEntriesVoxel) return;
1967 0 : TVectorF zres(7);
1968 0 : voxRes.flags = 0;
1969 0 : if (!TStatToolkit::LTMUnbinned(np,dz,zres,fLTMCut)) return;
1970 : //
1971 0 : float ab[2],err[3];
1972 0 : float sigMAD = FitPoly1Robust(np,tg,dy,ab,err,fLTMCut);
1973 0 : if (sigMAD<0) return;
1974 0 : float corrErr = err[0]*err[2];
1975 0 : corrErr = corrErr>0 ? err[1]/TMath::Sqrt(corrErr) : -999;
1976 : //printf("N:%3d A:%+e B:%+e / %+e %+e %+e | %+e %+e / %+e %+e\n",np,a,b,err[0],err[1],err[2], zres[1],zres[2], zres[3],zres[4]);
1977 : //
1978 0 : voxRes.D[kResX] = -ab[1];
1979 0 : voxRes.D[kResY] = ab[0];
1980 0 : voxRes.D[kResZ] = zres[1];
1981 0 : voxRes.E[kResX] = TMath::Sqrt(err[2]);
1982 0 : voxRes.E[kResY] = TMath::Sqrt(err[0]);
1983 0 : voxRes.E[kResZ] = zres[4];
1984 0 : voxRes.EXYCorr = corrErr;
1985 0 : voxRes.D[kResD] = voxRes.dYSigMAD = sigMAD; // later will be overriden by real dispersion
1986 0 : voxRes.dZSigLTM = zres[2];
1987 : //
1988 : // store the statistics
1989 0 : ULong64_t binStat = GetBin2Fill(voxRes.bvox,kVoxV);
1990 0 : voxRes.stat[kVoxV] = fArrNDStat[voxRes.bsec]->At(binStat);
1991 0 : for (int iv=kVoxDim;iv--;) voxRes.stat[iv] = fArrNDStat[voxRes.bsec]->At(binStat+iv-kVoxV);
1992 : //
1993 0 : voxRes.flags |= kDistDone;
1994 0 : }
1995 :
1996 : //_________________________________________________
1997 : void AliTPCDcalibRes::ProcessVoxelDispersions(int np, const float* tg, float *dy, bres_t& voxRes)
1998 : {
1999 : // extract Y (Z ignored at the moment) dispersions of the voxel
2000 : // correct Y distortions
2001 0 : if (np<2) return;
2002 0 : for (int i=np;i--;) dy[i] -= voxRes.DS[kResY] - voxRes.DS[kResX]*tg[i];
2003 0 : voxRes.D[kResD] = MAD2Sigma(np,dy);
2004 0 : voxRes.E[kResD] = voxRes.D[kResD]/TMath::Sqrt(2.*np); // a la gaussian RMS error, this is very crude
2005 0 : voxRes.flags |= kDispDone;
2006 : //
2007 0 : }
2008 :
2009 : //_____________________________________________________________________
2010 : Double_t AliTPCDcalibRes::GetLogL(TH1F* histo, int bin0, int bin1, double &mu, double &sig, double &logL0)
2011 : {
2012 : // Calculate log likelihood of normal distribution for the histo between boundaries
2013 : // bin0 and bin for given mu and sigma assumption. Exact Poisson statistics is assumed
2014 : // Also the approximate "reference" log-likelihood logLO is calculated in the following way:
2015 : // if the Poisson prob. for given bin is "m", then the logL0 gets contribution for this bin
2016 : // log("reference probability"), which I define as a geometric mean of probabilities for
2017 : // observing [M] and [M]+1 entries if m is large:
2018 : // P_ref = exp(-m)/m! M^m sqrt( m/(M+1) )
2019 : // -> ln(P_ref) = -(1/2+M)*ln(M/m) + M-m + 1/2 ln(M+1) - 1/2 ln(2 pi) -> -1/2 ln(2 pi m) for m>>1 (take m>5)
2020 : // (precise up to 1/2 ln(2 pi m) term of Stirling formula
2021 : // or = -m + m*log(m) - log( Gamma(1.+m) ) for m<~1
2022 : //
2023 : // integral
2024 : const double kNuLarge = 5.0, kMinSig2BinH = 0.01;
2025 0 : double dxh = 0.5*histo->GetBinWidth(1);
2026 0 : if ((sig/dxh)<kMinSig2BinH) {
2027 0 : AliWarningClassF("Too small sigma %.4e is provided for bin width %.4e",sig,dxh);
2028 0 : logL0 = -1;
2029 0 : return -1e9;
2030 : }
2031 : double sum=0, sum1=0, sum2=0;
2032 :
2033 0 : for (int ib=bin0;ib<=bin1;ib++) {
2034 0 : double w = histo->GetBinContent(ib);
2035 0 : double x = histo->GetBinCenter(ib);
2036 0 : sum += w;
2037 0 : sum1 += w*x;
2038 0 : sum2 += w*x*x;
2039 : }
2040 : //
2041 0 : double xb0 = histo->GetBinCenter(bin0)-dxh;
2042 0 : double xb1 = histo->GetBinCenter(bin1)+dxh;
2043 : //
2044 0 : if (sum<1e-6) {logL0 = -1e6; return -1e9;}
2045 0 : mu = sum1/sum;
2046 0 : sig = sum2/sum - mu*mu;
2047 0 : sig = TMath::Max(sig>0 ? TMath::Sqrt(sig) : 0.f, dxh/TMath::Sqrt(3)); // don't allow too small sigma
2048 :
2049 : //printf("Sample mu : %e sig: %e in %e %e\n",mu,sig,xb0,xb1);
2050 :
2051 : // estimated sig, mu are from the truncated sample, try to recover the truth
2052 0 : GetTruncNormMuSig(xb0,xb1, mu, sig);
2053 : //
2054 0 : xb0 -= mu;
2055 0 : xb1 -= mu;
2056 0 : double sqri2 = 1./(TMath::Sqrt(2.)*sig);
2057 : // normalization constant
2058 0 : double norm = 2.*sum / (TMath::Erf(xb1*sqri2) - TMath::Erf(xb0*sqri2));
2059 : //
2060 : // printf("Norm: %e\n",norm);
2061 : // likelihood
2062 : double logL = 0;
2063 0 : logL0 = 0;
2064 : const double kMinExp = 1e-100;
2065 0 : for (int i=bin0;i<=bin1;i++) {
2066 0 : double x = histo->GetBinCenter(i)-mu;
2067 0 : double w = histo->GetBinContent(i);
2068 0 : xb0 = x-dxh;
2069 0 : xb1 = x+dxh;
2070 : // bin expectation: normal integral within the bin
2071 0 : double nu = 0.5*norm*(TMath::Erf(xb1*sqri2) - TMath::Erf(xb0*sqri2));
2072 0 : if (nu<kMinExp) nu = kMinExp;
2073 0 : double logNFac = w<100 ? TMath::Log(TMath::Factorial(w)) : w*TMath::Log(w)-w + TMath::Log( sqrt(2*TMath::Pi()*w));
2074 0 : double logNu = TMath::Log(nu);
2075 0 : double logc = -nu + w*logNu - logNFac; // contribution of this bin to log-likelihood
2076 0 : logL += logc;
2077 : // now get the reference contribution
2078 : double logc0 = 0;
2079 0 : if (nu>kNuLarge) logc0 = -0.5*TMath::Log(2.*TMath::Pi()*nu);
2080 : else {
2081 0 : logc0 = -nu + nu*logNu - TMath::Log( TMath::Gamma(1.+nu) );
2082 : }
2083 0 : logL0 += logc0; // reference LL update
2084 : //printf("b: %d x:%+.2e nstd:%+.2e Exp:%e Obs:%e logc: %e logc0: %e\n",i,x,(x-mu)/sig, nu,w,logc, logc0);
2085 :
2086 : }
2087 : // printf("LogL: %e LogL0: %e\n",logL,logL0);
2088 : //
2089 : return logL;
2090 0 : }
2091 :
2092 : //___________________________________________________________________________
2093 : void AliTPCDcalibRes::TruncNormMod(double a, double b, double mu0, double sig0, double &muCf, double &sigCf)
2094 : {
2095 : // calculate truncated mean and sigma of normal distribution as
2096 : // mu_tr = mu0 + sig0*muCf
2097 : // sig_tr = sig0 * sigCf
2098 : //
2099 0 : const double sqrt2PiI = 1./TMath::Sqrt(TMath::Pi()*2.);
2100 0 : double sigI = 1./(sig0*TMath::Sqrt(2.));
2101 0 : double ra = (a-mu0)*sigI, rb = (b-mu0)*sigI;
2102 0 : double ra2 = ra*ra, rb2 = rb*rb;
2103 0 : double af = ra2<100 ? sqrt2PiI*TMath::Exp(-ra2) : 0;
2104 0 : double bf = rb2<100 ? sqrt2PiI*TMath::Exp(-rb2) : 0;
2105 : // double aF = 0.5*(1.+TMath::Erf(ra)), bF = 0.5*(1.+TMath::Erf(rb)), deltaF = bF-aF
2106 0 : double deltaF = 0.5*( TMath::Erf(rb) - TMath::Erf(ra) );
2107 0 : double deltaf = af - bf;
2108 0 : muCf = deltaf / deltaF;
2109 0 : sigCf = 1./TMath::Sqrt(1. + TMath::Sqrt(2)*(ra*af-rb*bf)/deltaF - muCf*muCf);
2110 : //
2111 0 : }
2112 :
2113 : //_____________________________________________________
2114 : Bool_t AliTPCDcalibRes::GetTruncNormMuSig(double a, double b, double &mean, double &sig)
2115 : {
2116 : // get estimate of real mu and sigma of normal distribution provided
2117 : // the mean and rms of sample truncated between a and b
2118 : const double kMinWindow=1e-2,kEpsRMS = 1e-4, kEpsMu = 1e-4;
2119 : const int kMaxIter = 200;
2120 : //
2121 0 : if (sig<1e-12) {
2122 0 : AliWarningClassF("Input sigma %e is too small",sig);
2123 0 : return kFALSE;
2124 : }
2125 0 : if ( (b-a)/sig<kMinWindow ) {
2126 0 : AliWarningClassF("Truncation window %e-%e is %e sigma only",a,b,(b-a)/sig);
2127 0 : return kFALSE;
2128 : }
2129 : //
2130 0 : double sig0=sig, mean0=mean; // initial values
2131 : // for protection, don't allow the sigma to grow above a factor of the flat distribution
2132 0 : double sigMax = 1.2*(b-a)/TMath::Sqrt(12.);
2133 : //
2134 0 : double m = mean, s = sig;
2135 0 : for (int i=0;i<kMaxIter;i++) {
2136 0 : double sclRMS,sclMU;
2137 0 : TruncNormMod(a,b,m,s, sclMU,sclRMS);
2138 : //
2139 0 : s = sig * sclRMS;
2140 : double mPrev = m, sPrev = s;
2141 0 : m = mean - sclMU*s;
2142 : //printf("%d -> M: %e S: %e\n",i,m, s);
2143 0 : if ( s>sigMax) {
2144 :
2145 : // printf("Iteration took sigma to twice of the flat distribution for "
2146 : // "mu0=%+.3e sig0=%.3e in %+.3e:%+.3e interval\n",mean0,sig0, a,b);
2147 0 : if (TMath::Abs(m-mean0)>sig0) {
2148 : // printf("Abandoning and returning input sigma and mean\n");
2149 : m = mean0; s = sig0;
2150 0 : }
2151 0 : break;
2152 : }
2153 0 : if (TMath::Abs(1.-sPrev/s)<kEpsRMS && TMath::Abs(m-mPrev)<kEpsMu ) break;
2154 0 : }
2155 : //
2156 0 : mean = m;
2157 0 : sig = s;
2158 :
2159 : return kTRUE;
2160 0 : }
2161 :
2162 : //_________________________________________________
2163 : void AliTPCDcalibRes::InitFieldGeom(Bool_t field,Bool_t geom)
2164 : {
2165 : // init geometry and field
2166 : // this requires the field and the geometry ...
2167 0 : if (fRun<1) AliFatal("Run number is not provided");
2168 0 : AliCDBManager* man = AliCDBManager::Instance();
2169 0 : if (fOCDBPath.IsNull()) fOCDBPath = "raw://";
2170 0 : if (!man->IsDefaultStorageSet()) man->SetDefaultStorage(fOCDBPath);
2171 0 : if (man->GetRun()!=fRun) man->SetRun(fRun);
2172 : //
2173 0 : Bool_t geomOK = AliGeomManager::GetGeometry() != 0;
2174 0 : if (geom && !geomOK) {
2175 0 : AliGeomManager::LoadGeometry();
2176 0 : AliGeomManager::ApplyAlignObjsFromCDB("TPC");
2177 0 : }
2178 0 : AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2179 0 : if (!fld && field) {
2180 0 : AliGRPManager grpMan;
2181 0 : grpMan.ReadGRPEntry();
2182 0 : grpMan.SetMagField();
2183 0 : fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2184 0 : }
2185 0 : if (fld) fBz = fld->SolenoidField();
2186 0 : }
2187 :
2188 :
2189 : //___________________________________________________________________________
2190 : THnF* AliTPCDcalibRes::CreateVoxelStatHisto(int sect)
2191 : {
2192 : // prepare histogram to store the means of distributions within the voxel
2193 :
2194 : // create binning for voxels statistics histograms
2195 0 : Int_t voxNBins[kVoxHDim];
2196 0 : Double_t voxBinMin[kVoxHDim],voxBinMax[kVoxHDim];
2197 0 : TString voxAxisName[kVoxHDim];
2198 :
2199 0 : voxAxisName[kVoxF] = "Y2X_Bin";
2200 0 : voxNBins[kVoxF] = fNY2XBins;
2201 0 : voxBinMin[kVoxF] = 0;
2202 0 : voxBinMax[kVoxF] = fNY2XBins;
2203 : //
2204 0 : voxAxisName[kVoxX] = "X_Bin";
2205 0 : voxNBins[kVoxX] = fNXBins;
2206 0 : voxBinMin[kVoxX] = 0;
2207 0 : voxBinMax[kVoxX] = fNXBins;
2208 : //
2209 0 : voxAxisName[kVoxZ] = "Z2X_Bin";
2210 0 : voxNBins[kVoxZ] = fNZ2XBins;
2211 0 : voxBinMin[kVoxZ] = 0;
2212 0 : voxBinMax[kVoxZ] = fNZ2XBins;
2213 : //
2214 0 : voxAxisName[kVoxV] = "Stat_Bin";
2215 0 : voxNBins[kVoxV] = kVoxHDim;
2216 0 : voxBinMin[kVoxV] = 0;
2217 0 : voxBinMax[kVoxV] = kVoxHDim;
2218 : //
2219 0 : THnF* h = new THnF(Form("hs%d",sect),"",kVoxHDim,voxNBins,voxBinMin,voxBinMax);
2220 0 : for (int i=0;i<kVoxHDim;i++) h->GetAxis(i)->SetName(voxAxisName[i].Data());
2221 0 : h->SetEntries(1); // otherwise drawing does not work well
2222 : return h;
2223 0 : }
2224 :
2225 : //===============================================================
2226 : //
2227 : // TRACK VALIDATION >>>>>>>>>>>>>>>>>>>>>>>>>>>
2228 : //
2229 :
2230 : //__________________________________________________________________________________
2231 : Bool_t AliTPCDcalibRes::ValidateTrack()
2232 : {
2233 : // if (nCl<fMinNCl) return kFALSE;
2234 0 : if (fNCl<fNVoisinMALong) return kFALSE;
2235 :
2236 0 : Bool_t rejCl[kNPadRows];
2237 0 : float rmsLong = 0.f;
2238 0 : int nRej = CheckResiduals(rejCl, rmsLong);
2239 0 : if (float(nRej)/fNCl > fMaxRejFrac) return kFALSE;
2240 0 : if (rmsLong>fMaxRMSLong) return kFALSE;
2241 : //
2242 : // flag outliers
2243 0 : for (int i=fNCl;i--;) if (rejCl[i]) fArrR[i] = fArrX[i] = -1;
2244 :
2245 0 : return kTRUE;
2246 0 : }
2247 :
2248 : //______________________________________________
2249 : void AliTPCDcalibRes::FixAlignmentBug(int sect, float q2pt, float bz, float& alp,
2250 : float& x, float &z, float &deltaY, float &deltaZ)
2251 : {
2252 : // fix alignment bug: https://alice.its.cern.ch/jira/browse/ATO-339?focusedCommentId=170850&page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel#comment-170850
2253 : //
2254 : // NOTE: deltaZ in the buggy code is calculated as Ztrack_with_bug - Zcluster_w/o_bug
2255 : static TGeoHMatrix *mCache[72] = {0};
2256 0 : if (sect<0||sect>=72) {
2257 0 : AliErrorF("Invalid sector %d",sect);
2258 0 : return;
2259 : }
2260 0 : int lr = sect/36 ? (AliGeomManager::kTPC2) : (AliGeomManager::kTPC1);
2261 0 : TGeoHMatrix* mgt = mCache[sect];
2262 0 : if (!mgt) {
2263 0 : int volID = AliGeomManager::LayerToVolUIDSafe(lr,sect%36);
2264 0 : mgt = new TGeoHMatrix(*AliGeomManager::GetTracking2LocalMatrix(volID));
2265 0 : mgt->MultiplyLeft(AliGeomManager::GetMatrix(volID));
2266 0 : mCache[sect] = mgt;
2267 0 : AliInfoF("Caching matrix for sector %d",sect);
2268 0 : }
2269 0 : double alpSect = ((sect%18)+0.5)*20.*TMath::DegToRad();
2270 :
2271 : // cluster in its proper alpha frame with alignment bug
2272 0 : double xyzClUse[3] = {x,0,z}; // this is what we read from the residual tree
2273 0 : double xyzTrUse[3] = {x, deltaY, z}; // track in bad cluster frame
2274 : //
2275 : // recover cluster Z position by adding deltaZ
2276 0 : double zClSave = xyzClUse[2] -= deltaZ; // here the cluster is not affected by Z alignment component of the bug!
2277 0 : static AliExternalTrackParam trDummy;
2278 0 : trDummy.Local2GlobalPosition(xyzClUse,alp); // misaligned cluster in global frame
2279 0 : double xyz0[3]={xyzClUse[0],xyzClUse[1],xyzClUse[2]};
2280 0 : mgt->MasterToLocal(xyz0,xyzClUse);
2281 : // we got ideal cluster in the sector tracking frame, but now the Z is wrong, since it was not affected by the bug!!!
2282 : //
2283 0 : xyzClUse[2] = zClSave;
2284 :
2285 : // go to ideal cluster frame
2286 0 : trDummy.Local2GlobalPosition(xyzClUse,alpSect); // ideal global
2287 0 : double alpFix = TMath::ATan2(xyzClUse[1],xyzClUse[0]); // fixed cluster phi
2288 0 : trDummy.Global2LocalPosition(xyzClUse,alpFix); // fixed cluster in in its frame
2289 : //
2290 0 : trDummy.Local2GlobalPosition(xyzTrUse,alp); // track in global frame
2291 0 : trDummy.Global2LocalPosition(xyzTrUse,alpFix); // track in cluster frame
2292 0 : alp = alpFix;
2293 : //
2294 0 : double dx = xyzTrUse[0] - xyzClUse[0]; // x might not be the same after alignment fix
2295 : // deduce track slopes assuming it comes from the vertex
2296 0 : double tgphi = tgpXY(xyzClUse[0],xyzTrUse[1]/xyzClUse[0],q2pt,bz);
2297 0 : xyzTrUse[1] -= dx*tgphi;
2298 0 : xyzTrUse[2] -= dx*xyzClUse[2]/xyzClUse[0]; // z2x
2299 : //
2300 0 : x = xyzClUse[0];
2301 0 : z = xyzTrUse[2]; // we still use track Z as a reference ...
2302 0 : deltaY = xyzTrUse[1]-xyzClUse[1];
2303 0 : deltaZ = xyzTrUse[2]-xyzClUse[2];
2304 : //
2305 0 : }
2306 :
2307 :
2308 : //_______________________________________________________________
2309 : int AliTPCDcalibRes::CheckResiduals(Bool_t* mask, float &rmsLongMA)
2310 : {
2311 :
2312 : int ip0=0,ip1;
2313 0 : int sec0 = fArrSectID[ip0];
2314 0 : int npLast = fNCl-1;
2315 : //
2316 : const int nMinAcc = 30;
2317 0 : float yDiffLL[kNPadRows] = {0.f};
2318 0 : float zDiffLL[kNPadRows] = {0.f};
2319 0 : float absDevY[kNPadRows] = {0.f};
2320 0 : float absDevZ[kNPadRows] = {0.f};
2321 :
2322 0 : rmsLongMA = 0.f;
2323 :
2324 0 : memset(mask,0,fNCl*sizeof(Bool_t));
2325 0 : for (int i=0;i<fNCl;i++) {
2326 0 : if (fArrSectID[i]==sec0 && i<npLast) continue;
2327 : //
2328 : // sector change or end of input reached
2329 : // run estimators for the points in the same sector
2330 0 : int npSec = i-ip0;
2331 0 : if (i==npLast) npSec++;
2332 : //
2333 0 : DiffToLocLine(npSec, fArrX+ip0, fArrDY+ip0, fNVoisinMA, yDiffLL+ip0);
2334 0 : DiffToLocLine(npSec, fArrX+ip0, fArrDZ+ip0, fNVoisinMA, zDiffLL+ip0);
2335 : // DiffToMA(npSec, fArrX+ip0, fArrDY+ip0, fNVoisinMA, yDiffLL+ip0);
2336 : // DiffToMA(npSec, fArrX+ip0, fArrDZ+ip0, fNVoisinMA, zDiffLL+ip0);
2337 : //
2338 : ip0 = i;
2339 0 : sec0 = fArrSectID[ip0];
2340 0 : }
2341 : // store abs deviations
2342 : int naccY=0,naccZ=0;
2343 0 : for (int i=fNCl;i--;) {
2344 0 : if (yDiffLL[i]) absDevY[naccY++] = TMath::Abs(yDiffLL[i]);
2345 0 : if (zDiffLL[i]) absDevZ[naccZ++] = TMath::Abs(zDiffLL[i]);
2346 : }
2347 : //
2348 : // estimate rms on 90% smallest deviations
2349 0 : int kmnY = 0.9*naccY,kmnZ = 0.9*naccZ;
2350 0 : if (naccY<nMinAcc || naccZ<nMinAcc) { // mask all
2351 0 : for (int i=fNCl;i--;) mask[i] = kTRUE;
2352 0 : return fNCl;
2353 : }
2354 :
2355 0 : SelKthMin(kmnY, naccY, absDevY);
2356 0 : SelKthMin(kmnZ, naccZ, absDevZ);
2357 : float rmsKY=0,rmsKZ=0;
2358 0 : for (int i=kmnY;i--;) rmsKY += absDevY[i]*absDevY[i];
2359 0 : for (int i=kmnZ;i--;) rmsKZ += absDevZ[i]*absDevZ[i];
2360 0 : rmsKY = TMath::Sqrt(rmsKY/kmnY);
2361 0 : rmsKZ = TMath::Sqrt(rmsKZ/kmnZ);
2362 : //
2363 0 : if (rmsKY<1e-6 || rmsKZ<1e-6) {
2364 0 : AliWarningF("Too small RMS: %f %f",rmsKY,rmsKZ);
2365 0 : for (int i=fNCl;i--;) mask[i] = kTRUE;
2366 0 : return fNCl;
2367 : }
2368 : //
2369 : // printf("RMSY %d min of %d: %f | RMSZ %d min of %d: %f\n",kmnY,naccY,rmsKY, kmnZ,naccZ,rmsKZ);
2370 : //
2371 : //
2372 0 : float rmsKYI = 1./rmsKY;
2373 0 : float rmsKZI = 1./rmsKZ;
2374 : int nMask=0, nacc = 0;
2375 0 : float yacc[kNPadRows],yDiffLong[kNPadRows];
2376 0 : for (int ip=0;ip<fNCl;ip++) {
2377 :
2378 0 : yDiffLL[ip] *= rmsKYI;
2379 0 : zDiffLL[ip] *= rmsKZI;
2380 0 : float dy = yDiffLL[ip], dz = zDiffLL[ip];
2381 0 : if (dy*dy+dz*dz>fMaxStdDevMA) {
2382 0 : mask[ip] = kTRUE;
2383 0 : nMask++;
2384 0 : }
2385 0 : else yacc[nacc++] = fArrDY[ip];
2386 : }
2387 : // rms to long-range moving average wrt surviving clusters
2388 0 : if (nacc>fNVoisinMALong) {
2389 0 : DiffToMA(nacc, yacc, fNVoisinMALong, yDiffLong);
2390 : float av=0,rms=0;
2391 0 : for (int i=0;i<nacc;i++) {
2392 0 : av += yDiffLong[i];
2393 0 : rms += yDiffLong[i]*yDiffLong[i];
2394 : }
2395 0 : av /= nacc;
2396 0 : rmsLongMA = rms/nacc - av*av;
2397 0 : rmsLongMA = rmsLongMA>0 ? TMath::Sqrt(rmsLongMA) : 0.f;
2398 0 : }
2399 : return nMask;
2400 : //
2401 0 : }
2402 :
2403 : //____________________________________________________________________
2404 : void AliTPCDcalibRes::FitCircle(int np, const float* x, const float* y,
2405 : double &xc, double &yc, double &r, float* dy)
2406 : {
2407 : // fit points to circle, if dy!=0, fill residuals
2408 : double x0=0.,y0=0.;
2409 0 : for (int i=np;i--;) {
2410 0 : x0 += x[i];
2411 0 : y0 += y[i];
2412 : }
2413 0 : x0 /= np;
2414 0 : y0 /= np;
2415 : double su2=0,sv2=0,suv=0,su3=0,sv3=0,su2v=0,suv2=0;
2416 0 : for (int i=np;i--;) {
2417 0 : double ui = x[i]-x0, ui2 = ui*ui;
2418 0 : double vi = y[i]-y0, vi2 = vi*vi;
2419 0 : suv += ui*vi;
2420 0 : su2 += ui2;
2421 0 : sv2 += vi2;
2422 0 : su3 += ui2*ui;
2423 0 : sv3 += vi2*vi;
2424 0 : su2v += ui2*vi;
2425 0 : suv2 += ui*vi2;
2426 : }
2427 0 : double rhsU = 0.5*(su3+suv2), rhsV = 0.5*(sv3+su2v);
2428 0 : double det = su2*sv2-suv*suv;
2429 0 : double uc = (rhsU*sv2 - rhsV*suv)/det;
2430 0 : double vc = (su2*rhsV - suv*rhsU)/det;
2431 0 : double r2 = uc*uc + vc*vc + (su2+sv2)/np;
2432 0 : xc = uc + x0;
2433 0 : yc = vc + y0;
2434 : //
2435 0 : if (dy) {
2436 0 : for (int i=np;i--;) {
2437 0 : double dx = x[i]-xc;
2438 0 : double dxr = r2 - dx*dx;
2439 0 : double ys = dxr>0 ? TMath::Sqrt(dxr) : 0;
2440 0 : double dy0 = y[i]-yc;
2441 0 : double dysp = dy0-ys;
2442 0 : double dysm = dy0+ys;
2443 0 : dy[i] = TMath::Abs(dysp)<TMath::Abs(dysm) ? dysp : dysm;
2444 : }
2445 0 : }
2446 0 : r = TMath::Sqrt(r2);
2447 0 : }
2448 :
2449 : //_________________________________________
2450 : void AliTPCDcalibRes::DiffToMA(int np, const float *y, const int winLR, float* diffMA)
2451 : {
2452 : // difference to moving average, excluding central element
2453 : //
2454 0 : double arrSumO[kNPadRows+1], *arrSum=arrSumO+1;
2455 0 : arrSum[-1] = 0.;
2456 0 : for (int ip=0;ip<np;ip++) arrSum[ip] = arrSum[ip-1]+y[ip];
2457 0 : for (int ip=0;ip<np;ip++) {
2458 0 : diffMA[ip] = 0;
2459 0 : int ipmn = ip-winLR;
2460 0 : int ipmx = ip+winLR;
2461 0 : if (ipmn<0) ipmn=0;
2462 0 : if (ipmx>=np) ipmx=np-1;
2463 0 : int nrm = (ipmx-ipmn);
2464 0 : if (nrm<winLR) continue;
2465 0 : double ma = (arrSum[ipmx]-arrSum[ipmn-1] - (arrSum[ip]-arrSum[ip-1]))/nrm;
2466 0 : diffMA[ip] = y[ip] - ma;
2467 0 : }
2468 :
2469 0 : }
2470 :
2471 :
2472 : //_______________________________________________
2473 : int AliTPCDcalibRes::DiffToLocLine(int np, const float* x, const float *y, const int nVoisin, float *diffY)
2474 : {
2475 : // calculate the difference between the point and linear extrapolation from neigbourhood
2476 : const float kEps = 1e-9;
2477 0 : double sumX1b[kNPadRows+1],sumX2b[kNPadRows+1],sumYb[kNPadRows+1],sumYXb[kNPadRows+1];
2478 0 : double *sumX1 = sumX1b+1, *sumX2 = sumX2b+1, *sumY0 = sumYb+1, *sumYX = sumYXb+1;
2479 : //
2480 0 : sumX1[-1]=0.f;
2481 0 : sumX2[-1]=0.f;
2482 0 : sumY0[-1]=0.f;
2483 0 : sumYX[-1]=0.f;
2484 : // accumulate matrix elements for whole array, element -1 is at 0
2485 : double sX1=0., sX2=0., sY0=0., sYX=0.;
2486 :
2487 0 : for (int ip=0;ip<np;ip++) {
2488 0 : sumX1[ip] = sX1 += x[ip];
2489 0 : sumX2[ip] = sX2 += x[ip]*x[ip];
2490 0 : sumY0[ip] = sY0 += y[ip];
2491 0 : sumYX[ip] = sYX += y[ip]*x[ip];
2492 0 : diffY[ip] = 0.0f;
2493 : }
2494 : //
2495 : int nAcc = 0;
2496 0 : for (int ip=0;ip<np;ip++) {
2497 0 : float &yEst = diffY[ip];
2498 : // estimate from the left
2499 0 : int ip0 = ip-nVoisin;
2500 0 : int ip1 = ip+nVoisin;
2501 0 : if (ip0<0) ip0=0;
2502 0 : if (ip1>=np) ip1=np-1;
2503 0 : int nrm = (ip1-ip0);
2504 0 : if (nrm<nVoisin) continue;
2505 0 : int ip0m = ip0-1, ipm = ip-1;
2506 : // extract sum from ip0 to ip1 from cumulant, S00=nrm, excluding tested point
2507 0 : sX1 = sumX1[ip1] - sumX1[ip0m] - (sumX1[ip]-sumX1[ipm]); // S01
2508 0 : sX2 = sumX2[ip1] - sumX2[ip0m] - (sumX2[ip]-sumX2[ipm]); // S11
2509 0 : sY0 = sumY0[ip1] - sumY0[ip0m] - (sumY0[ip]-sumY0[ipm]); // RHS0
2510 0 : sYX = sumYX[ip1] - sumYX[ip0m] - (sumYX[ip]-sumYX[ipm]); // RHS1
2511 0 : double det = nrm*sX2 - sX1*sX1;
2512 0 : if (det<kEps) continue;
2513 0 : double detI = 1./det;
2514 : // yLEst = offs + slop*x[ip]
2515 : // with offs=(sY0*sX2-sYX*sX1)/det and slop=(nVoisin*sYX-sY0*sX1)/det
2516 : // inverse err^2 = 1/(errOffs+x^2*errSlop+2*x*errSlpOffs) with
2517 : // errOffs = S11/det, errSlp=S00/det, errSlpOffs=-S01/det
2518 0 : yEst = y[ip]-((sY0*sX2-sYX*sX1) + (nrm*sYX-sY0*sX1)*x[ip])*detI;
2519 0 : nAcc++;
2520 : //
2521 0 : }
2522 0 : return nAcc;
2523 : //
2524 0 : }
2525 :
2526 : //====================================================================
2527 : int AliTPCDcalibRes::DiffToMedLine(int np, const float* x, const float *y, const int nVoisin, float *diffY)
2528 : {
2529 : int nAcc = 0;
2530 0 : float offs=0,slp=0;
2531 : float buff[kNPadRows];
2532 0 : for (int ip=0;ip<np;ip++) {
2533 0 : float &yEst = diffY[ip];
2534 0 : yEst = 0;
2535 : // estimate from the left
2536 0 : int ip0 = ip-nVoisin;
2537 0 : int ip1 = ip+nVoisin;
2538 0 : if (ip0<0) ip0=0;
2539 0 : if (ip1>=np) ip1=np-1;
2540 0 : int nrm = (ip1-ip0+1);
2541 0 : if (nrm<nVoisin) continue;
2542 : int ip0m = ip0-1, ipm = ip-1;
2543 0 : const float *arrx = x+ip0, *arry = y+ip0;
2544 0 : medFit(nrm, arrx, arry, offs,slp);
2545 : /*
2546 : float asum = 0;
2547 : for (int i=nrm;i--;) {
2548 : buff[i] = arry[i] - (offs + slp*arrx[i]);
2549 : if (i+ip0!=ip) asum += TMath::Abs(buff[i]);
2550 : }
2551 : asum /= nrm-1;
2552 : yEst = buff[ip-ip0]/asum;
2553 : */
2554 0 : yEst = y[ip] - (offs + slp*x[ip]);
2555 0 : nAcc++;
2556 0 : }
2557 0 : return nAcc;
2558 0 : }
2559 :
2560 : //_________________________________________________________
2561 : Int_t* AliTPCDcalibRes::LTMUnbinnedSig(int np, const float *arr, TVectorF ¶ms , Float_t sigTgt, Float_t minFrac, Bool_t sorted)
2562 : {
2563 : //
2564 : // LTM : Trimmed keeping at most minFrac of unbinned array to reach targer sigma
2565 : //
2566 : // Robust statistic to estimate properties of the distribution
2567 : // To handle binning error special treatment
2568 : // for definition of unbinned data see:
2569 : // http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999
2570 : //
2571 : // Function parameters:
2572 : // np - number of points in the array
2573 : // arr - data array (unsorted)
2574 : // params - vector with parameters
2575 : // - 0 - area
2576 : // - 1 - mean
2577 : // - 2 - rms
2578 : // - 3 - error estimate of mean
2579 : // - 4 - error estimate of RMS
2580 : // - 5 - first accepted element (of sorted array)
2581 : // - 6 - last accepted element (of sorted array)
2582 : //
2583 : // On success returns index of sorted events
2584 : //
2585 : static int *index = 0, book = 0;
2586 : static double* w = 0;
2587 0 : params[0] = 0.0f;
2588 0 : if (book<np) {
2589 0 : delete[] index;
2590 0 : book = np;
2591 0 : index = new int[book];
2592 0 : delete[] w;
2593 0 : w = new double[book+book];
2594 0 : }
2595 : //
2596 0 : double *wx1 = w, *wx2 = wx1+np;
2597 0 : if (!sorted) TMath::Sort(np,arr,index,kFALSE); // sort in increasing order
2598 0 : else for (int i=0;i<np;i++) index[i]=i;
2599 : // build cumulants
2600 : double sum1=0.0,sum2=0.0;
2601 0 : for (int i=0;i<np;i++) {
2602 0 : double x = arr[index[i]];
2603 0 : wx1[i] = (sum1+=x);
2604 0 : wx2[i] = (sum2+=x*x);
2605 : }
2606 : //
2607 : int keepMax = np;
2608 0 : int keepMin = minFrac*np;
2609 0 : if (keepMin>keepMax) keepMin = keepMax;
2610 : //
2611 0 : float sig2Tgt = sigTgt*sigTgt;
2612 0 : while (1) {
2613 0 : double minRMS = sum2+1e6;
2614 0 : int keepN = (keepMax+keepMin)>>1;
2615 0 : if (keepN<2) return 0;
2616 : //
2617 0 : params[0] = keepN;
2618 0 : int limI = np - keepN+1;
2619 0 : for (int i=0;i<limI;i++) {
2620 0 : int limJ = i+keepN-1;
2621 0 : Double_t sum1 = wx1[limJ] - (i ? wx1[i-1] : 0.0);
2622 0 : Double_t sum2 = wx2[limJ] - (i ? wx2[i-1] : 0.0);
2623 0 : double mean = sum1/keepN;
2624 0 : double rms2 = sum2/keepN - mean*mean;
2625 0 : if (rms2>minRMS) continue;
2626 : minRMS = rms2;
2627 0 : params[1] = mean;
2628 0 : params[2] = rms2;
2629 0 : params[5] = i;
2630 0 : params[6] = limJ;
2631 0 : }
2632 0 : if (minRMS<sig2Tgt) keepMin = keepN;
2633 : else keepMax = keepN;
2634 0 : if (keepMin>=keepMax-1) break;
2635 0 : }
2636 : //
2637 0 : if (!params[0]) return 0;
2638 0 : params[2] = TMath::Sqrt(params[2]);
2639 0 : params[3] = params[2]/TMath::Sqrt(params[0]); // error on mean
2640 0 : params[4] = params[3]/TMath::Sqrt(2.0); // error on RMS
2641 0 : return index;
2642 0 : }
2643 :
2644 : //___________________________________________________________________
2645 : float AliTPCDcalibRes::MAD2Sigma(int np, float* y)
2646 : {
2647 : // Sigma calculated from median absolute deviations, https://en.wikipedia.org/wiki/Median_absolute_deviation
2648 : // the input array is not modified
2649 0 : if (np<2) return 0;
2650 0 : int nph = np>>1;
2651 0 : if (nph&0x1) nph -= 1;
2652 : // don't abuse stack
2653 0 : float *ycHeap=0, ycStack[np<kMaxOnStack ? np:1],*yc=np<kMaxOnStack ? &ycStack[0] : (ycHeap = new float[np]);
2654 0 : memcpy(yc,y,np*sizeof(float));
2655 0 : float median = (np&0x1) ? SelKthMin(nph,np,yc) : 0.5f*(SelKthMin(nph-1,np,yc)+SelKthMin(nph,np,yc));
2656 : // build abs differences to median
2657 0 : for (int i=np;i--;) yc[i] = TMath::Abs(yc[i]-median);
2658 : // now get median of abs deviations
2659 0 : median = (np&0x1) ? SelKthMin(nph,np,yc) : 0.5f*(SelKthMin(nph-1,np,yc)+SelKthMin(nph,np,yc));
2660 0 : delete[] ycHeap; // if any...
2661 0 : return median*1.4826; // convert to Gaussian sigma
2662 0 : }
2663 :
2664 : //___________________________________________________________________
2665 : void AliTPCDcalibRes::medFit(int np, const float* x, const float* y, float &a, float &b, float* err,float delI)
2666 : {
2667 : // Median linear fit: minimizes abs residuals instead of squared ones
2668 : // Adapted from "Numerical Recipes in C"
2669 0 : float aa,bb,b1,b2,f,f1,f2,sigb,chisq=0.0f;
2670 0 : if (np<2) {
2671 0 : a = b = 0.0;
2672 0 : if (err) {
2673 0 : err[0] = err[1] = err[2] = 999.;
2674 0 : }
2675 0 : return;
2676 : }
2677 0 : if (!delI) {
2678 : float sx=0.0f,sxx=0.0f,sy=0.0f,sxy=0.0f,del;
2679 : //
2680 0 : for (int j=np;j--;) { sx += x[j]; sxx += x[j]*x[j];}
2681 0 : del = np*sxx-sx*sx;
2682 : //
2683 0 : for (int j=np;j--;) { sy += y[j]; sxy += x[j]*y[j];}
2684 : //
2685 0 : delI = 1./del;
2686 0 : aa = (sxx*sy-sx*sxy)*delI;
2687 0 : bb = (np*sxy-sx*sy)*delI;
2688 0 : if (err) {
2689 0 : err[0] = sxx*delI;
2690 0 : err[1] = sx*delI;
2691 0 : err[2] = np*delI;
2692 0 : }
2693 0 : }
2694 : else { // initial values provided
2695 0 : aa = a;
2696 0 : bb = b;
2697 : }
2698 : //
2699 0 : for (int j=np;j--;) {
2700 0 : float temp = y[j]-(aa+bb*x[j]);
2701 0 : chisq += temp*temp;
2702 : }
2703 : //
2704 0 : sigb = TMath::Sqrt(chisq*delI);
2705 : b1=bb;
2706 0 : f1 = RoFunc(np,x,y,b1,aa);
2707 0 : if (sigb>0) {
2708 0 : b2 = bb+TMath::Sign(float(3.0f*sigb),f1);
2709 0 : f2 = RoFunc(np,x,y,b2,aa);
2710 0 : if (f1==f2) {
2711 0 : a = aa;
2712 0 : b = bb;
2713 0 : return;
2714 : }
2715 0 : while (f1*f2 > 0.0f) { // bracketing
2716 0 : bb = b2 + 1.6f*(b2-b1);
2717 : b1 = b2;
2718 : f1 = f2;
2719 : b2 = bb;
2720 0 : f2 = RoFunc(np,x,y,b2,aa);
2721 : }
2722 0 : sigb = 0.01*sigb;
2723 0 : while (fabs(b2-b1)>sigb) {
2724 0 : bb = b1 + 0.5f*(b2-b1);
2725 0 : if (bb==b1 || bb==b2) break;
2726 0 : f = RoFunc(np,x,y,bb,aa);
2727 0 : if (f*f1 >= 0.0f) {
2728 : f1=f;
2729 : b1=bb;
2730 0 : }
2731 : else {
2732 : f2 = f;
2733 : b2 = bb;
2734 : }
2735 : }
2736 : }
2737 0 : a = aa;
2738 0 : b = bb;
2739 : //
2740 0 : }
2741 :
2742 : //___________________________________________________________________
2743 : float AliTPCDcalibRes::RoFunc(int np, const float* x, const float* y, float b, float &aa)
2744 : {
2745 : const float kEPS = 1.0e-7f;
2746 : static float* arrTmp = 0;
2747 : static int nBook = 0;
2748 0 : if (np>nBook) { // make sure the buffer is ok
2749 0 : nBook = np;
2750 0 : delete[] arrTmp;
2751 0 : arrTmp = new float[nBook];
2752 0 : }
2753 : float d,sum=0.0f;
2754 0 : for (int j=np;j--;) arrTmp[j] = y[j]-b*x[j];
2755 : //
2756 0 : int nph = np>>1;
2757 0 : if (np<20) { // it is faster to do insertion sort
2758 0 : for (int i=1;i<np;i++) {
2759 0 : float v = arrTmp[i];
2760 : int j;
2761 0 : for (j=i;j--;) if (arrTmp[j]>v) arrTmp[j+1]=arrTmp[j]; else break;
2762 0 : arrTmp[j+1] = v;
2763 : }
2764 0 : aa = (np&0x1) ? arrTmp[nph] : 0.5f*(arrTmp[nph-1]+arrTmp[nph]);
2765 0 : }
2766 : else {
2767 0 : aa = (np&0x1) ? SelKthMin(nph,np,arrTmp) :
2768 0 : 0.5f*(SelKthMin(nph-1,np,arrTmp)+SelKthMin(nph,np,arrTmp));
2769 : }
2770 0 : for (int j=np;j--;) {
2771 0 : d = y[j] - (b*x[j] + aa);
2772 0 : if (y[j] != 0.0f) d /= TMath::Abs(y[j]);
2773 0 : if (TMath::Abs(d) > kEPS) sum += (d >= 0.0f ? x[j] : -x[j]);
2774 : }
2775 0 : return sum;
2776 : }
2777 :
2778 : //___________________________________________________________________
2779 : Float_t AliTPCDcalibRes::SelKthMin(int k, int np, float* arr)
2780 : {
2781 : // Returns the k th smallest value in the array. The input array will be rearranged
2782 : // to have this value in location arr[k] , with all smaller elements moved before it
2783 : // (in arbitrary order) and all larger elements after (also in arbitrary order).
2784 : // From Numerical Recipes in C++
2785 :
2786 : int i,ir,j,l,mid;
2787 : float a;
2788 0 : l=0;ir=np-1;
2789 0 : for (;;) {
2790 0 : if (ir<=l+1) {
2791 0 : if (ir==l+1 && arr[ir]<arr[l]) swap(arr[l],arr[ir]);
2792 0 : return arr[k];
2793 : }
2794 : else {
2795 0 : int mid = (l+ir)>>1, i=l+1;
2796 0 : swap(arr[mid],arr[i]);
2797 0 : if (arr[i]>arr[ir]) swap(arr[i],arr[ir]);
2798 0 : if (arr[l]>arr[ir]) swap(arr[l] ,arr[ir]);
2799 0 : if (arr[i]>arr[l]) swap(arr[i],arr[l]);
2800 : j=ir;
2801 0 : a=arr[l];
2802 0 : for (;;) {
2803 0 : do i++; while (arr[i]<a);
2804 0 : do j--; while (arr[j]>a);
2805 0 : if (j<i) break;
2806 0 : swap(arr[i],arr[j]);
2807 : }
2808 0 : arr[l]=arr[j];
2809 0 : arr[j]=a;
2810 0 : if (j>=k) ir=j-1;
2811 0 : if (j<=k) l=i;
2812 : }
2813 : }
2814 : }
2815 :
2816 : //===============================================================
2817 :
2818 : //_________________________________________
2819 : void AliTPCDcalibRes::MakeVDriftOCDB(TString targetOCDBstorage)
2820 : {
2821 : // write OCDB object for VDrift (simplified copy/paste of AliTPCcalibAlignInterpolation::MakeVDriftOCDB)
2822 0 : TObjArray * driftArray = new TObjArray();
2823 : //
2824 0 : AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2825 0 : if (!fld) InitFieldGeom(kTRUE,kFALSE);
2826 0 : Double_t atime[2]={0,0};
2827 0 : Double_t deltaZ[2]={0,0};
2828 0 : Double_t t0[2]={0,0};
2829 0 : Double_t vdgy[2]={0,0};
2830 : //
2831 0 : atime[0] = fVDriftGraph->GetX()[0];
2832 0 : atime[1] = fVDriftGraph->GetX()[fVDriftGraph->GetN()-1];
2833 : //
2834 0 : AliTPCParam* param = AliTPCcalibDB::Instance()->GetParameters(); // just to get Vdrift
2835 : //
2836 0 : for (Int_t ipoint=0; ipoint<=1; ipoint++){
2837 0 : deltaZ[ipoint]=-(*fVDriftParam)[1]; // unit OK
2838 0 : vdgy[ipoint]=-(*fVDriftParam)[2]; // units OK
2839 0 : t0[ipoint]=-(*fVDriftParam)[0]/(1+(*fVDriftParam)[3]); // t0 to be normalized to the ms
2840 0 : t0[ipoint]/=(param->GetDriftV()/1000000.);
2841 : }
2842 0 : Double_t vdrifCorrRun=1./(1.+(*fVDriftParam)[3]);
2843 : //
2844 : TGraphErrors *graphDELTAZ=0;
2845 : TGraphErrors *graphT0=0;
2846 : TGraphErrors *graphVDGY=0;
2847 : TGraphErrors *graphDRIFTVD=0;
2848 : TGraphErrors *graph=0;
2849 : //
2850 0 : TString grPrefix="ALIGN_ITSB_TPC_";
2851 : //
2852 0 : graphDELTAZ=new TGraphErrors(2, atime, deltaZ);
2853 0 : graphDELTAZ->SetName(grPrefix+"DELTAZ");
2854 0 : driftArray->AddLast(graphDELTAZ);
2855 0 : graphT0=new TGraphErrors(2, atime, t0);
2856 0 : graphT0->SetName(grPrefix+"T0");
2857 0 : driftArray->AddLast(graphT0);
2858 0 : graphVDGY=new TGraphErrors(2, atime, vdgy);
2859 0 : graphVDGY->SetName(grPrefix+"VDGY");
2860 0 : driftArray->AddLast(graphVDGY);
2861 : //
2862 : // drift velocity
2863 0 : graph = new TGraphErrors(*fVDriftGraph);
2864 0 : graph->SetName(grPrefix+"DRIFTVD");
2865 0 : Int_t npoints = graph->GetN();
2866 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++) {
2867 0 : graph->GetY()[ipoint] = (vdrifCorrRun*(1-graph->GetY()[ipoint]))-1.;
2868 0 : graph->GetEY()[ipoint] = vdrifCorrRun*graph->GetEY()[ipoint];
2869 : }
2870 0 : driftArray->AddLast(graph);
2871 : //
2872 : AliCDBStorage* targetStorage = 0x0;
2873 0 : if (targetOCDBstorage.Length()==0) {
2874 0 : targetOCDBstorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
2875 0 : targetStorage = AliCDBManager::Instance()->GetStorage(targetOCDBstorage.Data());
2876 0 : }
2877 0 : else if (targetOCDBstorage.CompareTo("same",TString::kIgnoreCase) == 0 ){
2878 0 : targetStorage = AliCDBManager::Instance()->GetDefaultStorage();
2879 0 : }
2880 : else {
2881 0 : targetStorage = AliCDBManager::Instance()->GetStorage(targetOCDBstorage.Data());
2882 : }
2883 :
2884 0 : AliCDBMetaData* metaData = new AliCDBMetaData;
2885 0 : metaData->SetObjectClassName("TObjArray");
2886 0 : metaData->SetResponsible("Marian Ivanov");
2887 0 : metaData->SetBeamPeriod(1);
2888 0 : metaData->SetAliRootVersion(">v5-07-20"); //AliRoot version
2889 0 : metaData->SetComment("AliTPCcalibAlignInterpolation Calibration of the time dependence of the drift velocity using Residual trees");
2890 0 : AliCDBId id1("TPC/Calib/TimeDrift", fRun, fRun);
2891 :
2892 : //now the entry owns the metadata, but NOT the data
2893 0 : AliCDBEntry *driftCDBentry=new AliCDBEntry(driftArray,id1,metaData,kFALSE);
2894 0 : targetStorage->Put(driftCDBentry);
2895 : //
2896 0 : delete driftCDBentry;
2897 0 : }
2898 :
2899 : //_________________________________________
2900 : void AliTPCDcalibRes::LoadVDrift()
2901 : {
2902 : // load vdrift params
2903 0 : fVDriftGraph = 0;
2904 0 : fVDriftParam = 0;
2905 0 : TString vdname = Form("%s.root",kDriftFileName);
2906 : TFile *fdrift = 0;
2907 0 : if (gSystem->AccessPathName(vdname.Data()) || !(fdrift=TFile::Open(vdname.Data())) ) {
2908 0 : AliWarningF("vdrift file %s not accessible",vdname.Data());
2909 0 : return;
2910 : }
2911 0 : TTree * tree = (TTree*)fdrift->Get("fitTimeStat");
2912 0 : if (!tree) AliWarningF("tree fitTimeStat not avaliable in file %s.root",kDriftFileName);
2913 : else {
2914 0 : tree->SetBranchAddress("grTRDReg.",&fVDriftGraph);
2915 0 : tree->SetBranchAddress("paramRobust.",&fVDriftParam);
2916 0 : tree->GetEntry(0);
2917 0 : if (fVDriftGraph==NULL || fVDriftGraph->GetN()<=0) {
2918 0 : AliWarning("ITS/TRD drift calibration not availalble. Trying ITS/TOF");
2919 0 : tree->SetBranchAddress("grTOFReg.",&fVDriftGraph);
2920 0 : tree->GetEntry(0);
2921 : }
2922 : }
2923 0 : delete tree;
2924 0 : delete fdrift;
2925 0 : }
2926 :
2927 : //_________________________________________
2928 : float AliTPCDcalibRes::GetDriftCorrection(float z, float x, float phi, int rocID)
2929 : {
2930 : // apply vdrift correction
2931 0 : if (!fVDriftParam) return 0.f;
2932 0 : int side = ((rocID/kNSect)&0x1) ? -1:1; // C:A
2933 0 : float drift = side>0 ? kZLim[0]-z : z+kZLim[1];
2934 0 : float gy = TMath::Sin(phi)*x;
2935 : Double_t pvecFit[3];
2936 0 : pvecFit[0]= side; // z shift (cm)
2937 0 : pvecFit[1]= drift*gy/kZLim[side<0]; // global y gradient
2938 0 : pvecFit[2]= drift; // drift length
2939 0 : float expected = (fVDriftParam==NULL) ? 0:
2940 0 : (*fVDriftParam)[0]+
2941 0 : (*fVDriftParam)[1]*pvecFit[0]+
2942 0 : (*fVDriftParam)[2]*pvecFit[1]+
2943 0 : (*fVDriftParam)[3]*pvecFit[2];
2944 0 : return -side*(expected+fCorrTime*drift);
2945 : //
2946 0 : }
2947 :
2948 : //_____________________________________________________
2949 : float AliTPCDcalibRes::tgpXY(float x, float y, float q2p, float bz)
2950 : {
2951 : // get the tg of primary track inclination wrt padrow given
2952 : // that it was registered at X,Y sector coordinates
2953 0 : float c = q2p*bz*(-0.299792458e-3f);
2954 0 : if (TMath::Abs(c)<1e-9) return y/x;
2955 0 : float r2 = x*x+y*y;
2956 0 : float det = 4./r2 - c*c;
2957 : float snp = 0;
2958 0 : if (det<0) {
2959 0 : snp = TMath::Sign(-0.8f,c);
2960 0 : AliWarningF("track of q2p=%f cannot reach x:%f y:%f",q2p,x,y);
2961 0 : }
2962 : else {
2963 0 : snp = 0.5f*(y*TMath::Sqrt(det)-c*x); // snp at vertex
2964 0 : snp += x*c; // snp at x,y
2965 : }
2966 0 : return snp/TMath::Sqrt((1.f-snp)*(1.f+snp));
2967 0 : }
2968 :
2969 :
2970 : //=========================================================================
2971 : //
2972 : // IO related methods
2973 : //
2974 : //=========================================================================
2975 : //
2976 : //___________________________________________________________________
2977 : void AliTPCDcalibRes::WriteStatHistos()
2978 : {
2979 : // write stat histos
2980 0 : TString statOutName = Form("%s.root",kStatOut);
2981 0 : TFile* statOutFile = TFile::Open(statOutName.Data(),"recreate");
2982 0 : for (int is=0;is<kNSect2;is++) fStatHist[is]->Write("", TObject::kOverwrite);
2983 0 : statOutFile->Close();
2984 0 : delete statOutFile;
2985 : //
2986 0 : }
2987 :
2988 : //___________________________________________________________________
2989 : void AliTPCDcalibRes::LoadStatHistos()
2990 : {
2991 : // load bin stat histos
2992 0 : if (fStatHist[0]) return; // histos are in memory
2993 0 : TString statOutName = Form("%s.root",kStatOut);
2994 0 : TFile* statOutFile = TFile::Open(statOutName.Data());
2995 0 : if (!statOutFile) AliFatalF("LoadStatHistos: failed to read file %s",statOutName.Data());
2996 0 : for (int is=0;is<kNSect2;is++) {
2997 0 : fStatHist[is] = (THnF*) statOutFile->Get(Form("hs%d",is));
2998 0 : if (!fStatHist[is]) AliFatalF("LoadStatHistos: failed to read secto %d histo from %s",is,statOutName.Data());
2999 0 : fArrNDStat[is] = (TNDArrayT<float>*)&fStatHist[is]->GetArray();
3000 : }
3001 0 : statOutFile->Close();
3002 0 : delete statOutFile;
3003 : //
3004 0 : }
3005 :
3006 : //___________________________________________________________________
3007 : void AliTPCDcalibRes::WriteResTree()
3008 : {
3009 : // output file for results tree
3010 0 : TStopwatch sw;
3011 0 : sw.Start();
3012 0 : bres_t voxRes, *voxResP=&voxRes;
3013 :
3014 0 : AliSysInfo::AddStamp("ResTree",0,0,0,0);
3015 0 : if (fChebCorr) fChebCorr->Init();
3016 0 : TFile* flOut = new TFile(GetVoxResFileName(),"recreate");
3017 0 : TTree* resTree = new TTree("voxRes","final distortions, see GetListOfAliases");
3018 0 : resTree->Branch("res", &voxRes);
3019 0 : resTree->Branch("L",&fLumiCOG);
3020 0 : for (int i=0;i<kVoxDim;i++) {
3021 0 : resTree->SetAlias(kVoxName[i],Form("bvox[%d]",i));
3022 0 : resTree->SetAlias(Form("%sAV",kVoxName[i]),Form("stat[%d]",i));
3023 : }
3024 0 : resTree->SetAlias(Form("%sAV",kVoxName[kVoxV]),Form("stat[%d]",kVoxV));
3025 0 : for (int i=0;i<kResDim;i++) {
3026 0 : resTree->SetAlias(kResName[i],Form("D[%d]",i));
3027 0 : resTree->SetAlias(Form("%sS",kResName[i]),Form("DS[%d]",i));
3028 0 : resTree->SetAlias(Form("%sC",kResName[i]),Form("DC[%d]",i));
3029 0 : resTree->SetAlias(Form("%sE",kResName[i]),Form("E[%d]",i));
3030 : }
3031 : //
3032 0 : resTree->SetAlias("fitOK",Form("(flags&0x%x)==0x%x",kDistDone,kDistDone));
3033 0 : resTree->SetAlias("dispOK",Form("(flags&0x%x)==0x%x",kDispDone,kDispDone));
3034 0 : resTree->SetAlias("smtOK",Form("(flags&0x%x)==0x%x",kSmoothDone,kSmoothDone));
3035 0 : resTree->SetAlias("masked",Form("(flags&0x%x)==0x%x",kMasked,kMasked));
3036 : //
3037 0 : for (int is=0;is<kNSect2;is++) {
3038 :
3039 0 : bres_t* sectData = fSectGVoxRes[is];
3040 0 : if (!sectData) {
3041 0 : AliWarningF("No processed data for sector %d",is);
3042 0 : continue;
3043 : }
3044 :
3045 : // now store the data for the sector
3046 0 : for (int iz=0;iz<fNZ2XBins;iz++) {
3047 0 : for (int ix=0;ix<fNXBins;ix++) {
3048 0 : for (int ip=0;ip<fNY2XBins;ip++) {
3049 0 : int binGlo = GetVoxGBin(ix,ip,iz);
3050 0 : bres_t *voxel = §Data[binGlo];
3051 0 : if (fChebCorr) {
3052 0 : int row = GetRowID(voxel->stat[kVoxX]);
3053 : int roc = is;
3054 0 : if (row>=kNRowIROC) {roc += kNSect2; row -= kNRowIROC;}
3055 0 : fChebCorr->Eval(roc,row, voxel->stat[kVoxF],voxel->stat[kVoxZ],voxel->DC);
3056 0 : }
3057 0 : memcpy(&voxRes,voxel,sizeof(bres_t)); // store in the sector data array
3058 0 : resTree->Fill();
3059 : }
3060 : }
3061 : }
3062 : //
3063 0 : } // end of sector loop
3064 : //
3065 : //
3066 : // write small trees with luminosity and vdrift info
3067 : //
3068 0 : TTree* infoTree = new TTree("runInfo","run information");
3069 0 : if (fLumiCTPGraph) infoTree->Branch("lumiCTP", &fLumiCTPGraph);
3070 0 : if (fVDriftGraph) infoTree->Branch("driftGraph", &fVDriftGraph);
3071 0 : if (fVDriftParam) infoTree->Branch("driftParam", &fVDriftParam);
3072 0 : infoTree->Branch("tminCTPrun",&fTMinCTP);
3073 0 : infoTree->Branch("tmaxCTPrun",&fTMaxCTP);
3074 0 : infoTree->Branch("tminTBin",&fTMin);
3075 0 : infoTree->Branch("tmaxTBin",&fTMax);
3076 0 : infoTree->Fill();
3077 : //
3078 0 : flOut->cd();
3079 0 : resTree->Write("", TObject::kOverwrite);
3080 0 : infoTree->Write("", TObject::kOverwrite);
3081 0 : delete resTree;
3082 0 : delete infoTree;
3083 : //
3084 0 : if (fChebCorr) {
3085 0 : fChebCorr->Write();
3086 : }
3087 : //
3088 0 : flOut->Close();
3089 0 : delete flOut;
3090 : //
3091 0 : sw.Stop();
3092 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
3093 0 : AliSysInfo::AddStamp("ResTree",1,0,0,0);
3094 :
3095 0 : }
3096 :
3097 : //___________________________________________________________________
3098 : Bool_t AliTPCDcalibRes::LoadDataFromResTree(const char* resTreeFile)
3099 : {
3100 : // Fill voxels info from existing resVox tree for reprocessing
3101 0 : TStopwatch sw;
3102 0 : sw.Start();
3103 0 : bres_t voxRes, *voxResP=&voxRes;
3104 : //
3105 0 : TTree* resTree = LoadTree(resTreeFile);
3106 0 : if (!resTree) {AliErrorF("Failed to extract resTree from %s",resTreeFile); return kFALSE;}
3107 0 : resTree->SetBranchAddress("res", &voxResP);
3108 : //
3109 0 : int nent = resTree->GetEntries();
3110 0 : if (nent != kNSect2*fNGVoxPerSector) AliFatalF("resTree from %s has %d voxels per sector, this object: %d",
3111 : resTreeFile,nent/kNSect2,fNGVoxPerSector);
3112 : //
3113 0 : for (int is=0;is<kNSect2;is++) {
3114 0 : if (fSectGVoxRes[is]) delete[] fSectGVoxRes[is];
3115 0 : bres_t* sectData = fSectGVoxRes[is] = new bres_t[fNGVoxPerSector];
3116 : }
3117 0 : for (int ient=0;ient<nent;ient++) {
3118 : //
3119 0 : resTree->GetEntry(ient);
3120 0 : bres_t* sectData = fSectGVoxRes[voxRes.bsec];
3121 0 : int binGlo = GetVoxGBin(voxRes.bvox);
3122 0 : bres_t *voxel = §Data[binGlo];
3123 0 : memcpy(voxel,&voxRes,sizeof(bres_t));
3124 : // the X distortion contribution was already subtracted from the Y,Z components, restore it
3125 0 : voxel->D[kResZ] -= voxel->stat[kVoxZ]*voxel->DS[kResX];
3126 : // reset smoothed params
3127 0 : voxel->flags &= ~(kSmoothDone|kMasked);
3128 0 : for (int ir=kResDim;ir--;) voxel->DS[ir]=voxel->DC[ir]=0.f;
3129 : } // end of sector loop
3130 : //
3131 0 : CloseTreeFile(resTree);
3132 : //
3133 0 : sw.Stop();
3134 0 : AliInfoF("timing: real: %.3f cpu: %.3f",sw.RealTime(), sw.CpuTime());
3135 : //
3136 : return kTRUE;
3137 0 : }
3138 :
3139 :
3140 : //=========================================================================
3141 : //
3142 : // fitting/smoothing related methods
3143 : //
3144 : //=========================================================================
3145 :
3146 : //_____________________________________________________
3147 : Bool_t AliTPCDcalibRes::FitPoly2(const float* x,const float* y, const float* w, int np, float *res, float *err)
3148 : {
3149 : // poly2 fitter
3150 0 : if (np<3) return kFALSE; // no enough points
3151 : double sumW[5]={0},sumY[3]={0};
3152 0 : for (int ip=np;ip--;) {
3153 0 : double ww = w ? w[ip] : 1.0;
3154 0 : sumW[0] += ww;
3155 0 : sumY[0] += ww*y[ip];
3156 0 : sumW[1] += (ww*=x[ip]);
3157 0 : sumY[1] += ww*y[ip];
3158 0 : sumW[2] += (ww*=x[ip]);
3159 0 : sumY[2] += ww*y[ip];
3160 0 : sumW[3] += (ww*=x[ip]);
3161 0 : sumW[4] += (ww*=x[ip]);
3162 : }
3163 0 : double min00 = (sumW[2]*sumW[4]-sumW[3]*sumW[3]);
3164 0 : double min01 = (sumW[1]*sumW[4]-sumW[3]*sumW[2]);
3165 0 : double min02 = (sumW[1]*sumW[3]-sumW[2]*sumW[2]);
3166 0 : double min11 = (sumW[0]*sumW[4]-sumW[2]*sumW[2]);
3167 0 : double min12 = (sumW[0]*sumW[3]-sumW[2]*sumW[1]);
3168 0 : double min22 = (sumW[0]*sumW[2]-sumW[1]*sumW[1]);
3169 0 : double det = sumW[0]*min00
3170 0 : - sumW[1]*min01
3171 0 : + sumW[2]*min02;
3172 0 : if (TMath::Abs(det)<1e-12) return kFALSE;
3173 0 : double detI = 1./det;
3174 0 : double det0 = sumY[0]*min00
3175 0 : - sumW[1]*(sumY[1]*sumW[4]-sumW[3]*sumY[2])
3176 0 : + sumW[2]*(sumY[1]*sumW[3]-sumW[2]*sumY[2]);
3177 0 : double det1 = sumW[0]*(sumY[1]*sumW[4]-sumW[3]*sumY[2])
3178 0 : - sumY[0]*min01
3179 0 : + sumW[2]*(sumW[1]*sumY[2]-sumY[1]*sumW[2]);
3180 0 : double det2 = sumW[0]*(sumW[2]*sumY[2]-sumY[1]*sumW[3])
3181 0 : - sumW[1]*(sumW[1]*sumY[2]-sumY[1]*sumW[2])
3182 0 : + sumY[0]*min02;
3183 0 : res[0] = det0*detI;
3184 0 : res[1] = det1*detI;
3185 0 : res[2] = det2*detI;
3186 : //
3187 0 : if (err) {
3188 0 : err[0] = min00*detI; // e00
3189 0 : err[1] =-min01*detI; // e10
3190 0 : err[2] = min11*detI; // e11
3191 0 : err[3] = min02*detI; // e20
3192 0 : err[4] =-min12*detI; // e21
3193 0 : err[5] = min22*detI; // e21
3194 0 : }
3195 : //
3196 : return kTRUE;
3197 0 : }
3198 :
3199 : //_____________________________________________________
3200 : Bool_t AliTPCDcalibRes::FitPoly1(const float* x,const float* y, const float* w, int np, float *res, float *err)
3201 : {
3202 : // poly1 fitter
3203 0 : if (np<2) return kFALSE; // no enough points
3204 : double sumW[3]={0},sumY[2]={0};
3205 0 : for (int ip=np;ip--;) {
3206 0 : double ww = w ? w[ip]:1.0;
3207 0 : sumW[0] += ww;
3208 0 : sumY[0] += ww*y[ip];
3209 0 : sumW[1] += (ww*=x[ip]);
3210 0 : sumY[1] += ww*y[ip];
3211 0 : sumW[2] += (ww*=x[ip]);
3212 : }
3213 0 : double det = sumW[0]*sumW[2] - sumW[1]*sumW[1];
3214 0 : if (TMath::Abs(det)<1e-12) return kFALSE;
3215 0 : double detI = 1./det;
3216 0 : double det0 = sumY[0]*sumW[2] - sumW[1]*sumY[1];
3217 0 : double det1 = sumW[0]*sumY[1] - sumY[0]*sumW[1];
3218 0 : res[0] = det0*detI;
3219 0 : res[1] = det1*detI;
3220 : //
3221 0 : if (err) {
3222 0 : err[0] = sumW[2]*detI; // e00
3223 0 : err[1] =-sumW[1]*detI; // e10
3224 0 : err[2] = sumW[0]*detI; // e11
3225 0 : }
3226 : //
3227 : return kTRUE;
3228 0 : }
3229 :
3230 : //________________________________
3231 : Int_t AliTPCDcalibRes::ValidateVoxels(int isect)
3232 : {
3233 : // apply voxel validation cuts, calculate number of low stat or masked voxels
3234 : int cntMasked=0, cntInvalid=0;
3235 : //
3236 0 : fXBinIgnore[isect].ResetAllBits();
3237 0 : bres_t* sectData = fSectGVoxRes[isect];
3238 0 : for (int ix=0;ix<fNXBins;ix++) {
3239 : int nvalidXBin = 0;
3240 0 : for (int ip=0;ip<fNY2XBins;ip++) {
3241 0 : for (int iz=0;iz<fNZ2XBins;iz++) { // extract line in z
3242 0 : int binGlo = GetVoxGBin(ix,ip,iz);
3243 0 : bres_t *voxRes = §Data[binGlo];
3244 0 : Bool_t voxOK = voxRes->flags&kDistDone && !(voxRes->flags&kMasked);
3245 0 : if (voxOK) {
3246 : // check fit errors
3247 0 : if (voxRes->E[kResY]*voxRes->E[kResY]>fMaxFitYErr2 ||
3248 0 : voxRes->E[kResX]*voxRes->E[kResX]>fMaxFitXErr2 ||
3249 0 : TMath::Abs(voxRes->EXYCorr)>fMaxFitXYCorr) voxOK = kFALSE;
3250 : // check raw distributions sigmas
3251 0 : if (voxRes->dYSigMAD > fMaxSigY || voxRes->dZSigLTM > fMaxSigZ) voxOK = kFALSE;
3252 0 : if (!voxOK) cntMasked++;
3253 : }
3254 0 : if (voxOK) {
3255 0 : nvalidXBin++;
3256 0 : voxRes->flags |= kDistDone;
3257 0 : }
3258 : else {
3259 0 : cntInvalid++;
3260 0 : voxRes->flags |= kMasked;
3261 : }
3262 : } // loop over Z
3263 : } // loop over Y
3264 : //
3265 0 : fValidFracXBin[isect][ix] = Float_t(nvalidXBin)/(fNY2XBins*fNZ2XBins);
3266 0 : if (fValidFracXBin[isect][ix]<fMinValidVoxFracDrift) {
3267 0 : AliWarningF("Sector%2d: Xbin%3d has %4.1f%% of voxels valid (%d out of %d)",
3268 : isect,ix,100*fValidFracXBin[isect][ix],nvalidXBin,fNY2XBins*fNZ2XBins);
3269 0 : }
3270 : //
3271 : } // loop over X
3272 : //
3273 : // mask X-bins which cannot be smoothed
3274 : int nValidInPatch = 0;
3275 : // 1st loop: find bad regions
3276 0 : short nbadReg=0,badStart[kNPadRows],badEnd[kNPadRows];
3277 : Bool_t prevBad = kFALSE;
3278 : float fracBadRows = 0;
3279 0 : for (int ix=0;ix<fNXBins;ix++) {
3280 0 : if (fValidFracXBin[isect][ix]<fMinValidVoxFracDrift) {
3281 0 : fracBadRows++;
3282 0 : if (prevBad) badEnd[nbadReg] = ix;
3283 : else {
3284 0 : badStart[nbadReg] = badEnd[nbadReg] = ix;
3285 : prevBad = kTRUE;
3286 : }
3287 : }
3288 : else {
3289 0 : if (prevBad) {
3290 : prevBad = kFALSE;
3291 0 : nbadReg++;
3292 0 : }
3293 : }
3294 : }
3295 0 : if (prevBad) nbadReg++; // last bad region was not closed
3296 : //
3297 0 : fracBadRows /= fNXBins;
3298 0 : if (fracBadRows>fMaxBadRowsPerSector) {
3299 0 : AliWarningF("Sector%2d: Fraction of bad X-bins %.3f > %.3f: masking whole sector",
3300 : isect,fracBadRows,fMaxBadRowsPerSector);
3301 0 : for (int ix=0;ix<fNXBins;ix++) SetXBinIgnored(isect,ix);
3302 0 : }
3303 : else {
3304 : //
3305 : // 2nd loop: disable those regions which cannot be smoothed
3306 0 : for (int ibad=0;ibad<nbadReg;ibad++) {
3307 0 : short badSize=badEnd[ibad]-badStart[ibad]+1,badSizeNext=ibad<(nbadReg-1) ? badEnd[ibad]-badStart[ibad]+1 : 0;
3308 : // disable too large bad patches
3309 0 : if (badSize>fMaxBadXBinsToCover) for (int i=0;i<badSize;i++) SetXBinIgnored(isect,badStart[ibad]+i);
3310 : // disable too small isolated good patches
3311 0 : if (badSizeNext>fMaxBadXBinsToCover && (badStart[ibad+1]-badEnd[ibad]-1)<fMinGoodXBinsToCover) {
3312 0 : for (int i=badEnd[ibad]+1;i<badStart[ibad+1];i++) SetXBinIgnored(isect,i);
3313 0 : }
3314 : }
3315 0 : if (nbadReg) {
3316 : int ib=0;
3317 : // is 1st good patch too small?
3318 0 : if (GetXBinIgnored(isect,badStart[0]) && badStart[ib]<fMinGoodXBinsToCover) {
3319 0 : for (int i=0;i<badStart[ib];i++) SetXBinIgnored(isect,i);
3320 0 : }
3321 : // last good patch is too small?
3322 0 : ib = nbadReg-1;
3323 0 : if (GetXBinIgnored(isect,badStart[ib]) && (fNXBins-badEnd[ib]-1)<fMinGoodXBinsToCover) {
3324 0 : for (int i=badEnd[ib]+1;i<fNXBins;i++) SetXBinIgnored(isect,i);
3325 0 : }
3326 0 : }
3327 : //
3328 : }
3329 : //
3330 0 : int nMaskedRows = fXBinIgnore[isect].CountBits();
3331 0 : AliInfoF("Sector%2d: Voxel stat: Masked: %5d(%07.3f%%) Invalid: %5d(%07.3f%%) -> Masked %3d rows out of %3d",isect,
3332 : cntMasked, 100*float(cntMasked)/fNGVoxPerSector,
3333 : cntInvalid,100*float(cntInvalid)/fNGVoxPerSector,
3334 : nMaskedRows,fNXBins);
3335 : //
3336 0 : return fNXBins-nMaskedRows;
3337 0 : }
3338 :
3339 : //________________________________
3340 : Int_t AliTPCDcalibRes::Smooth0(int isect)
3341 : {
3342 : // apply linear regression kernel smoother
3343 : int cnt = 0;
3344 0 : bres_t* sectData = fSectGVoxRes[isect];
3345 0 : for (int ix=0;ix<fNXBins;ix++) {
3346 0 : if (GetXBinIgnored(isect,ix)) continue;
3347 0 : for (int ip=0;ip<fNY2XBins;ip++) {
3348 0 : for (int iz=0;iz<fNZ2XBins;iz++) { // extract line in z
3349 0 : int binGlo = GetVoxGBin(ix,ip,iz);
3350 0 : bres_t *vox = §Data[binGlo];
3351 0 : vox->flags &= ~kSmoothDone;
3352 0 : Bool_t res = GetSmoothEstimate(vox->bsec,vox->stat[kVoxX],vox->stat[kVoxF],vox->stat[kVoxZ],
3353 : BIT(kResX)|BIT(kResY)|BIT(kResZ), // at this moment we cannot smooth dispersion
3354 0 : vox->DS);
3355 0 : if (!res) fNSmoothingFailedBins[isect]++;
3356 0 : else vox->flags |= kSmoothDone;
3357 : }
3358 : }
3359 0 : }
3360 : // now subtract the dX contribution to DZ
3361 0 : for (int ix=0;ix<fNXBins;ix++) {
3362 0 : if (GetXBinIgnored(isect,ix)) continue;
3363 0 : for (int ip=0;ip<fNY2XBins;ip++) {
3364 0 : for (int iz=0;iz<fNZ2XBins;iz++) {
3365 0 : int binGlo = GetVoxGBin(ix,ip,iz);
3366 0 : bres_t *vox = §Data[binGlo];
3367 0 : if (!(vox->flags&kSmoothDone)) continue;
3368 0 : vox->DS[kResZ] += vox->stat[kVoxZ]*vox->DS[kResX]; // remove slope*dx contribution from DZ
3369 0 : vox->D[kResZ] += vox->stat[kVoxZ]*vox->DS[kResX]; // remove slope*dx contribution from DZ
3370 0 : }
3371 : }
3372 0 : }
3373 : //
3374 0 : return cnt;
3375 : }
3376 :
3377 : //________________________________________________________________
3378 : Bool_t AliTPCDcalibRes::GetSmoothEstimate(int isect, float x, float p, float z, Int_t which, float *res, float *deriv)
3379 : {
3380 : // get smooth estimate of distortions mentioned in "which" bit pattern for point in sector coordinates (x,y/x,z/x)
3381 : // smoothing results also saved in the fLastSmoothingRes (allow derivative calculation)
3382 : //
3383 0 : int minPointsDir[kVoxDim]={0}; // min number of points per direction
3384 : //
3385 : const float kTrialStep = 0.5;
3386 0 : Bool_t doDim[kResDim] = {kFALSE};
3387 0 : for (int i=0;i<kResDim;i++) {
3388 0 : doDim[i] = (which&(0x1<<i))>0;
3389 0 : if (doDim[i]) res[i] = 0;
3390 : }
3391 : //
3392 : // extimate smoothing matrix size and min number of points
3393 : int matSize = kSmtLinDim;
3394 0 : for (int i=0;i<kVoxDim;i++) {
3395 0 : minPointsDir[i] = 3; // for pol1 smoothing require at least 3 points
3396 0 : if (fSmoothPol2[i]) {
3397 0 : minPointsDir[i]++;
3398 0 : matSize++;
3399 0 : }
3400 : }
3401 0 : double cmat[kResDim][kMaxSmtDim*(kMaxSmtDim+1)/2];
3402 : static int maxNeighb = 10*10*10;
3403 0 : static bres_t **currClus = new bres_t*[maxNeighb];
3404 0 : static float* currCache = new float[maxNeighb*kVoxHDim];
3405 : //
3406 : //loop over neighbours which can contribute
3407 : //
3408 : //
3409 0 : int ix0,ip0,iz0;
3410 0 : FindVoxel(x,p, isect<kNSect ? z : -z, ix0,ip0,iz0); // find nearest voxel
3411 0 : bres_t* sectData = fSectGVoxRes[isect];
3412 0 : int binCen = GetVoxGBin(ix0,ip0,iz0); // global bin of nearest voxel
3413 0 : bres_t* voxCen = §Data[binCen]; // nearest voxel
3414 : //
3415 0 : int maxTrials[kVoxDim];
3416 0 : maxTrials[kVoxZ] = fNBins[kVoxZ]/2;
3417 0 : maxTrials[kVoxF] = fNBins[kVoxF]/2;
3418 0 : maxTrials[kVoxX] = fMaxBadXBinsToCover*2;
3419 :
3420 0 : int trial[kVoxDim]={0};
3421 0 : while(1) {
3422 : //
3423 0 : memset(fLastSmoothingRes,0,kResDim*kMaxSmtDim*sizeof(double));
3424 0 : memset(cmat,0,kResDim*kMaxSmtDim*(kMaxSmtDim+1)/2*sizeof(double));
3425 : //
3426 : int nbOK=0; // accounted neighbours
3427 : //
3428 0 : float stepX = fStepKern[kVoxX]*(1. + kTrialStep*trial[kVoxX]);
3429 0 : float stepF = fStepKern[kVoxF]*(1. + kTrialStep*trial[kVoxF]);
3430 0 : float stepZ = fStepKern[kVoxZ]*(1. + kTrialStep*trial[kVoxZ]);
3431 : //
3432 0 : if (!(voxCen->flags&kDistDone) || (voxCen->flags&kMasked) || GetXBinIgnored(isect,ix0)) {
3433 : // closest voxel has no data, increase smoothing step
3434 0 : stepX+=kTrialStep*fStepKern[kVoxX];
3435 0 : stepF+=kTrialStep*fStepKern[kVoxF];
3436 0 : stepZ+=kTrialStep*fStepKern[kVoxZ];
3437 0 : }
3438 : //
3439 : // effective kernel widths accounting for the increased bandwidth at the edges and missing data
3440 0 : float kWXI = GetDXI(ix0) *fKernelWInv[kVoxX]*fStepKern[kVoxX]/stepX;
3441 0 : float kWFI = GetDY2XI(ix0)*fKernelWInv[kVoxF]*fStepKern[kVoxF]/stepF;
3442 0 : float kWZI = GetDZ2XI() *fKernelWInv[kVoxZ]*fStepKern[kVoxZ]/stepZ;
3443 0 : int istepX = TMath::Nint(stepX+0.5);
3444 0 : int istepF = TMath::Nint(stepF+0.5);
3445 0 : int istepZ = TMath::Nint(stepZ+0.5);
3446 : // for edge bins increase kernel size and neighbours search
3447 0 : int ixMn=ix0-istepX,ixMx=ix0+istepX;
3448 0 : if (ixMn<0) {
3449 : ixMn = 0;
3450 0 : ixMx = TMath::Min(TMath::Nint(ix0+stepX*fKernelScaleEdge[kVoxX]),fNXBins-1);
3451 0 : kWXI /= fKernelScaleEdge[kVoxX];
3452 0 : }
3453 0 : if (ixMx>=fNXBins) {
3454 0 : ixMx = fNXBins-1;
3455 0 : ixMn = TMath::Max(TMath::Nint(ix0-stepX*fKernelScaleEdge[kVoxX]),0);
3456 0 : kWXI /= fKernelScaleEdge[kVoxX];
3457 0 : }
3458 : //
3459 0 : int ipMn=ip0-istepF,ipMx=ip0+istepF;
3460 0 : if (ipMn<0) {
3461 : ipMn = 0;
3462 0 : ipMx = TMath::Min(TMath::Nint(ip0+stepF*fKernelScaleEdge[kVoxF]),fNY2XBins-1);
3463 0 : kWFI /= fKernelScaleEdge[kVoxF];
3464 0 : }
3465 0 : if (ipMx>=fNY2XBins) {
3466 0 : ipMx = fNY2XBins-1;
3467 0 : ipMn = TMath::Max(TMath::Nint(ip0-stepF*fKernelScaleEdge[kVoxF]),0);
3468 0 : kWFI /= fKernelScaleEdge[kVoxF];
3469 0 : }
3470 : //
3471 0 : int izMn=iz0-istepZ,izMx=iz0+istepZ;
3472 0 : if (izMn<0) {
3473 : izMn = 0;
3474 0 : izMx = TMath::Min(TMath::Nint(iz0+stepZ*fKernelScaleEdge[kVoxZ]),fNZ2XBins-1);
3475 0 : kWZI /= fKernelScaleEdge[kVoxZ];
3476 0 : }
3477 0 : if (izMx>=fNZ2XBins) {
3478 0 : izMx = fNZ2XBins-1;
3479 0 : izMn = TMath::Max(TMath::Nint(iz0-stepZ*fKernelScaleEdge[kVoxZ]),0);
3480 0 : kWZI /= fKernelScaleEdge[kVoxZ];
3481 0 : }
3482 : //
3483 0 : int nOccX[ixMx-ixMn+1],nOccF[ipMx-ipMn+1],nOccZ[izMx-izMn+1];
3484 : //
3485 : // check if cache arrays should be expanded
3486 0 : int nbCheck = (ixMx-ixMn+1)*(ipMx-ipMn+1)*(izMx-izMn+1);
3487 0 : if (nbCheck>=maxNeighb) { // need to expand caches
3488 0 : int mxNb = nbCheck+100;
3489 0 : delete[] currClus;
3490 0 : delete[] currCache;
3491 0 : currClus = new bres_t*[mxNb];
3492 0 : currCache= new float[mxNb*kVoxHDim];
3493 0 : maxNeighb = mxNb;
3494 0 : }
3495 : //
3496 0 : for (int i=ixMx-ixMn+1;i--;) nOccX[i]=0;
3497 0 : for (int i=ipMx-ipMn+1;i--;) nOccF[i]=0;
3498 0 : for (int i=izMx-izMn+1;i--;) nOccZ[i]=0;
3499 0 : double u2Vec[3];
3500 : //1st loop, check presence of enough points, cache precalculated values
3501 0 : float *cacheVal = currCache;
3502 0 : for (int ix=ixMn;ix<=ixMx;ix++) {
3503 0 : for (int ip=ipMn;ip<=ipMx;ip++) {
3504 0 : for (int iz=izMn;iz<=izMx;iz++) {
3505 : //
3506 0 : int binNb = GetVoxGBin(ix,ip,iz); // global bin
3507 0 : bres_t* voxNb = §Data[binNb];
3508 0 : if (!(voxNb->flags&kDistDone) || (voxNb->flags&kMasked) || GetXBinIgnored(isect,ix)) continue; // skip voxels w/o data
3509 : // estimate weighted distance
3510 0 : float dx = voxNb->stat[kVoxX]-x;
3511 0 : float df = voxNb->stat[kVoxF]-p;
3512 0 : float dz = voxNb->stat[kVoxZ]-z;
3513 0 : float dxw = dx*kWXI, dfw = df*kWFI, dzw = dz*kWZI;
3514 0 : u2Vec[0] = dxw*dxw;
3515 0 : u2Vec[1] = dfw*dfw;
3516 0 : u2Vec[2] = dzw*dzw;
3517 0 : double kernW = GetKernelWeight(u2Vec,3);
3518 0 : if (kernW<kZeroK) continue;
3519 : // new point is validated
3520 0 : nOccX[ix-ixMn]++;
3521 0 : nOccF[ip-ipMn]++;
3522 0 : nOccZ[iz-izMn]++;
3523 0 : currClus[nbOK] = voxNb;
3524 0 : cacheVal[kVoxX] = dx;
3525 0 : cacheVal[kVoxF] = df;
3526 0 : cacheVal[kVoxZ] = dz;
3527 0 : cacheVal[kVoxV] = kernW;
3528 0 : cacheVal += kVoxHDim;
3529 0 : nbOK++;
3530 0 : }
3531 : }
3532 : }
3533 : //
3534 : // check if we have enough points in every dimension
3535 0 : int np[kVoxDim]={0};
3536 0 : for (int i=ixMx-ixMn+1;i--;) if (nOccX[i]) np[kVoxX]++;
3537 0 : for (int i=ipMx-ipMn+1;i--;) if (nOccF[i]) np[kVoxF]++;
3538 0 : for (int i=izMx-izMn+1;i--;) if (nOccZ[i]) np[kVoxZ]++;
3539 0 : Bool_t enoughPoints = kTRUE, incrDone[kVoxDim] = {0};
3540 0 : for (int i=0;i<kVoxDim;i++) {
3541 0 : if (np[i]<minPointsDir[i]) { // need to extend smoothing neighborhood
3542 : enoughPoints=kFALSE;
3543 0 : if (trial[i]<maxTrials[i] && !incrDone[i]) { //try to increment only missing direction
3544 0 : trial[i]++; incrDone[i]=kTRUE;
3545 0 : }
3546 0 : else if (trial[i]==maxTrials[i]) { // cannot increment missing direction, try others
3547 0 : for (int j=kVoxDim;j--;) {
3548 0 : if (i!=j && trial[j]<maxTrials[j] && !incrDone[j]) {
3549 0 : trial[j]++; incrDone[j]=kTRUE;
3550 0 : }
3551 : }
3552 0 : }
3553 : }
3554 : }
3555 0 : if (!enoughPoints) {
3556 0 : if (!(incrDone[kVoxX]||incrDone[kVoxF]||incrDone[kVoxZ])) {
3557 0 : AliErrorF("Voxel Z:%d F:%d X:%d Trials limit reached: Z:%d F:%d X:%d",
3558 : voxCen->bvox[kVoxZ],voxCen->bvox[kVoxF],voxCen->bvox[kVoxX],
3559 : trial[kVoxZ],trial[kVoxF],trial[kVoxX]);
3560 0 : return kFALSE;
3561 : }
3562 : /*
3563 : AliWarningF("Sector:%2d x=%.3f y/x=%.3f z/x=%.3f (iX:%d iY2X:%d iZ2X:%d)\n"
3564 : "not enough neighbours (need min %d) %d %d %d (tot: %d) | Steps: %.1f %.1f %.1f\n"
3565 : "trying to increase filter bandwidth (trialXFZ:%d %d %d)\n",
3566 : isect,x,p,z,ix0,ip0,iz0,2,np[kVoxX],np[kVoxF],np[kVoxZ],nbOK,stepX,stepF,stepZ,
3567 : trial[kVoxX],trial[kVoxF],trial[kVoxZ]);
3568 : */
3569 0 : continue;
3570 : }
3571 : //
3572 : // now fill the matrices and solve
3573 0 : cacheVal = currCache;
3574 0 : for (int ib=0;ib<nbOK;ib++) {
3575 0 : double kernW = cacheVal[kVoxV];
3576 0 : double dx=cacheVal[kVoxX], df=cacheVal[kVoxF], dz=cacheVal[kVoxZ],dx2=dx*dx, df2=df*df, dz2 = dz*dz;
3577 0 : cacheVal += kVoxHDim;
3578 0 : const bres_t* voxNb = currClus[ib];
3579 0 : for (int id=0;id<kResDim;id++) {
3580 0 : if (!doDim[id]) continue;
3581 : double kernWD = kernW;
3582 0 : if (fUseErrInSmoothing) kernWD /= (voxNb->E[id]*voxNb->E[id]); // apart from the kernel value, account for the point error
3583 0 : double *cmatD = cmat[id];
3584 0 : double *rhsD = &fLastSmoothingRes[id*kMaxSmtDim];
3585 : //
3586 0 : double kernWDx=kernWD*dx, kernWDf=kernWD*df, kernWDz=kernWD*dz;
3587 0 : double kernWDx2=kernWDx*dx, kernWDf2=kernWDf*df, kernWDz2=kernWDz*dz;
3588 : //
3589 : // linear part
3590 : int el=-1,elR=-1;
3591 0 : cmatD[++el] += kernWD;
3592 0 : rhsD[++elR] += kernWD*voxNb->D[id];
3593 : //
3594 0 : cmatD[++el] += kernWDx; cmatD[++el] += kernWDx2;
3595 0 : rhsD[++elR] += kernWDx*voxNb->D[id];
3596 : //
3597 0 : cmatD[++el] += kernWDf; cmatD[++el] += kernWDx*df; cmatD[++el] += kernWDf2;
3598 0 : rhsD[++elR] += kernWDf*voxNb->D[id];
3599 : //
3600 0 : cmatD[++el] += kernWDz; cmatD[++el] += kernWDx*dz; cmatD[++el] += kernWDf*dz; cmatD[++el] += kernWDz2;
3601 0 : rhsD[++elR] += kernWDz*voxNb->D[id];
3602 : //
3603 : // check if quadratic part is needed
3604 0 : if (fSmoothPol2[kVoxX]) {
3605 0 : cmatD[++el] += kernWDx2; cmatD[++el] += kernWDx2*dx; cmatD[++el] += kernWDf*dx2; cmatD[++el] += kernWDz*dx2; cmatD[++el] += kernWDx2*dx2;
3606 0 : rhsD[++elR] += kernWDx2*voxNb->D[id];
3607 0 : }
3608 0 : if (fSmoothPol2[kVoxF]) {
3609 0 : cmatD[++el] += kernWDf2; cmatD[++el] += kernWDx*df2; cmatD[++el] += kernWDf*df2; cmatD[++el] += kernWDz*df2; cmatD[++el] += kernWDx2*df2; cmatD[++el] += kernWDf2*df2;
3610 0 : rhsD[++elR] += kernWDf2*voxNb->D[id];
3611 0 : }
3612 0 : if (fSmoothPol2[kVoxZ]) {
3613 0 : cmatD[++el] += kernWDz2; cmatD[++el] += kernWDx*dz2; cmatD[++el] += kernWDf*dz2; cmatD[++el] += kernWDz*dz2; cmatD[++el] += kernWDx2*dz2; cmatD[++el] += kernWDf2*dz2; cmatD[++el] += kernWDz2*dz2;
3614 0 : rhsD[++elR] += kernWDz2*voxNb->D[id];
3615 0 : }
3616 0 : }
3617 : }
3618 : //
3619 : Bool_t fitRes = kTRUE;
3620 : //
3621 : // solve system of linear equations
3622 0 : AliSymMatrix mat(matSize);
3623 0 : for (int id=0;id<kResDim;id++) {
3624 0 : if (!doDim[id]) continue;
3625 0 : mat.Reset();
3626 0 : double *cmatD = cmat[id];
3627 0 : double *rhsD = &fLastSmoothingRes[id*kMaxSmtDim];
3628 : int el=-1,elR=-1,row=-1;
3629 0 : mat(++row,0) = cmatD[++el];
3630 0 : mat(++row,0) = cmatD[++el]; mat(row,1) = cmatD[++el];
3631 0 : mat(++row,0) = cmatD[++el]; mat(row,1) = cmatD[++el]; mat(row,2) = cmatD[++el];
3632 0 : mat(++row,0) = cmatD[++el]; mat(row,1) = cmatD[++el]; mat(row,2) = cmatD[++el]; mat(row,3) = cmatD[++el];
3633 : // pol2 elements if needed
3634 0 : if (fSmoothPol2[kVoxX]) {
3635 : const int colLim = (++row)+1;
3636 0 : for (int col=0;col<colLim;col++) mat(row,col) = cmatD[++el];
3637 0 : }
3638 0 : if (fSmoothPol2[kVoxF]) {
3639 0 : const int colLim = (++row)+1;
3640 0 : for (int col=0;col<colLim;col++) mat(row,col) = cmatD[++el];
3641 0 : }
3642 0 : if (fSmoothPol2[kVoxZ]) {
3643 0 : const int colLim = (++row)+1;
3644 0 : for (int col=0;col<colLim;col++) mat(row,col) = cmatD[++el];
3645 0 : }
3646 : //
3647 0 : fitRes &= mat.SolveChol(rhsD);
3648 0 : if (!fitRes) {
3649 0 : for (int i=kVoxDim;i--;) trial[i]++;
3650 0 : AliWarningF("Sector:%2d x=%.3f y/x=%.3f z/x=%.3f (iX:%d iY2X:%d iZ2X:%d)\n"
3651 : "neighbours range used %d %d %d (tot: %d) | Steps: %.1f %.1f %.1f\n"
3652 : "Solution for smoothing Failed, trying to increase filter bandwidth (trialXFZ: %d %d %d)",
3653 : isect,x,p,z,ix0,ip0,iz0,np[kVoxX],np[kVoxF],np[kVoxZ],nbOK,stepX,stepF,stepZ,trial[kVoxX],trial[kVoxF],trial[kVoxZ]);
3654 0 : continue;
3655 : }
3656 0 : res[id] = rhsD[0];
3657 0 : if (deriv) for (int j=0;j<3;j++) deriv[id*3 +j] = rhsD[j+1]; // ignore eventual pol2 term
3658 0 : }
3659 : //
3660 : break; // success
3661 0 : } // end of loop over allowed trials
3662 :
3663 0 : return kTRUE;
3664 :
3665 0 : }
3666 :
3667 : //_____________________________________________________________________
3668 : Bool_t AliTPCDcalibRes::GetSmooth1D
3669 : (double xQuery, // query X
3670 : double valerr[2], // output array for value and its error
3671 : int np, const double* x, const double* y, const double* err, // points x,y, optional error
3672 : double w, int kType,Bool_t usePol2, // kernel width and type, do we use pol1 or pol2 interpolation
3673 : Bool_t xIncreasing, // are points increasing in X
3674 : float maxD2Range // may increase kernel width up to maxD2Range*dataRange
3675 : ) const
3676 : {
3677 : // get kernel (width w) smoothed value at xQuery for (opionally ordered in x-increasing) array x,y (optionally err)
3678 0 : int minPoints = 3+usePol2;
3679 : // don't abuse stack
3680 0 : float *ySelHeap=0,ySelStack[np<kMaxOnStack ? np:1],*ySel=np<kMaxOnStack ? &ySelStack[0] : (ySelHeap=new float[np]);
3681 0 : float *xSelHeap=0,xSelStack[np<kMaxOnStack ? np:1],*xSel=np<kMaxOnStack ? &xSelStack[0] : (xSelHeap=new float[np]);
3682 0 : float *wSelHeap=0,wSelStack[np<kMaxOnStack ? np:1],*wSel=np<kMaxOnStack ? &wSelStack[0] : (wSelHeap=new float[np]);
3683 : //
3684 : // select only points which have chance to contribute
3685 : double ws = w;
3686 0 : double xmax = x[np-1], xmin = x[0];
3687 0 : if (!xIncreasing) {
3688 0 : for (int i=np;i--;) {
3689 0 : if (xmax>x[i]) xmax = x[i];
3690 0 : if (xmin<x[i]) xmin = x[i];
3691 : }
3692 0 : }
3693 0 : double maxD,w2Sum=0.,range = (xmax - xmin)*maxD2Range;
3694 : //
3695 : const float kEps = 1e-12;
3696 : int npUse=0;
3697 0 : do {
3698 0 : maxD = ws*(kType==kGaussianKernel ? kMaxGaussStdDev : 1.0f);
3699 0 : double ws2 = 1./(ws*ws);
3700 : npUse = 0;
3701 : w2Sum = 0.;
3702 0 : for (int ip=0;ip<np;ip++) {
3703 0 : double df = x[ip]-xQuery;
3704 0 : if (df>maxD) {
3705 0 : if (xIncreasing) break; // with sorted arrays no point in checking further
3706 0 : continue;
3707 : }
3708 0 : if (df<-maxD) continue;
3709 0 : xSel[npUse] = df; // we fit wrt queried point!
3710 0 : ySel[npUse] = y[ip];
3711 0 : df *= df*ws2;
3712 0 : wSel[npUse] = GetKernelWeight(&df,1);
3713 0 : w2Sum += wSel[npUse];
3714 0 : if (wSel[npUse]<kEps) continue;
3715 0 : if (err && err[ip]>0) wSel[npUse] /= err[ip]*err[ip];
3716 0 : npUse++;
3717 0 : }
3718 : //
3719 0 : if (npUse>=minPoints) break;
3720 0 : ws *= 1.5; // expand kernel
3721 0 : }
3722 0 : while (maxD<range);
3723 : //
3724 : Bool_t res = kFALSE;
3725 0 : if (npUse<minPoints) {
3726 0 : valerr[0] = valerr[1] = 0.f;
3727 0 : AliErrorF("Found only %d points, %d required",npUse,minPoints);
3728 0 : }
3729 : else {
3730 0 : float resVal[3],resErr[6];
3731 0 : if ((res=usePol2 ? FitPoly2(xSel,ySel,wSel,npUse,resVal,resErr) : FitPoly1(xSel,ySel,wSel,npUse,resVal,resErr))) {
3732 0 : valerr[0] = resVal[0];
3733 0 : valerr[1] = TMath::Sqrt(resErr[0]/w2Sum);
3734 0 : }
3735 0 : }
3736 0 : delete[] ySelHeap;
3737 0 : delete[] xSelHeap;
3738 0 : delete[] wSelHeap;
3739 : //
3740 0 : return res;
3741 0 : }
3742 :
3743 : //_____________________________________
3744 : Double_t AliTPCDcalibRes::GetKernelWeight(double* u2vec,int np) const
3745 : {
3746 : double w = 1;
3747 0 : if (fKernelType == kEpanechnikovKernel) {
3748 0 : for (int i=np;i--;) {
3749 0 : if (u2vec[i]>1) return 0.;
3750 0 : w *= 3./4.*(1.-u2vec[i]);
3751 : }
3752 : }
3753 0 : else if (fKernelType == kGaussianKernel) {
3754 : double u2 = 0;
3755 0 : for (int i=np;i--;) u2 += u2vec[i];
3756 0 : w = u2<kMaxGaussStdDev*kMaxGaussStdDev*np ? TMath::Exp(-u2)/TMath::Sqrt(2.*TMath::Pi()) : 0;
3757 0 : }
3758 : else {
3759 0 : AliFatalF("Kernel type %d is not defined",fKernelType);
3760 : }
3761 0 : return w;
3762 0 : }
3763 :
3764 :
3765 : //_____________________________________
3766 : void AliTPCDcalibRes::SetKernelType(int tp, float bwX, float bwP, float bwZ, float scX,float scP,float scZ)
3767 : {
3768 : // set kernel type and widths in terms of binning in X,Y/X and Z/X, define aux variables
3769 0 : fKernelType = tp;
3770 : //
3771 0 : fKernelScaleEdge[kVoxX] = scX;
3772 0 : fKernelScaleEdge[kVoxF] = scP;
3773 0 : fKernelScaleEdge[kVoxZ] = scZ;
3774 :
3775 0 : fKernelWInv[kVoxX] = bwX>0 ? 1./bwX : 1.;
3776 0 : fKernelWInv[kVoxF] = bwP>0 ? 1./bwP : 1.;
3777 0 : fKernelWInv[kVoxZ] = bwZ>0 ? 1./bwZ : 1.;
3778 :
3779 0 : if (fKernelType == kEpanechnikovKernel) { // bandwidth 1
3780 0 : fStepKern[kVoxX] = TMath::Nint(bwX+0.5);
3781 0 : fStepKern[kVoxF] = TMath::Nint(bwP+0.5);
3782 0 : fStepKern[kVoxZ] = TMath::Nint(bwZ+0.5);
3783 0 : }
3784 : else if (kGaussianKernel) { // look in ~5 sigma
3785 0 : fStepKern[kVoxX] = TMath::Nint(bwX*5.+0.5);
3786 0 : fStepKern[kVoxF] = TMath::Nint(bwP*5.+0.5);
3787 0 : fStepKern[kVoxZ] = TMath::Nint(bwZ*5.+0.5);
3788 : }
3789 : else {
3790 : AliFatalF("Kernel type %d is not defined",fKernelType);
3791 : }
3792 0 : for (int i=kVoxDim;i--;) if (fStepKern[i]<1) fStepKern[i] = 1;
3793 : //
3794 0 : }
3795 :
3796 :
3797 : //=========================================================================
3798 :
3799 : //_____________________________________________
3800 : void AliTPCDcalibRes::CreateCorrectionObject()
3801 : {
3802 : // create correction object for given time slice
3803 :
3804 0 : AliSysInfo::AddStamp("CreateCorrectionObject",0,0,0,0);
3805 : //
3806 : // check if there are failures
3807 : Bool_t create = kTRUE;
3808 0 : for (int i=0;i<kNSect2;i++) {
3809 0 : if (fNSmoothingFailedBins[i]>0) {
3810 0 : AliErrorF("%d failed voxels in sector %d",fNSmoothingFailedBins[i],i);
3811 : create = kFALSE;
3812 0 : }
3813 : }
3814 : //
3815 0 : if (!create) {
3816 0 : AliError("ATTENTION: MAP WILL NOT BE CREATED");
3817 0 : return;
3818 : }
3819 : //
3820 0 : TString name = Form("run%d_%lld_%lld",fRun,fTMin,fTMax);
3821 0 : fChebCorr = new AliTPCChebCorr(name.Data(),name.Data(),
3822 0 : fChebPhiSlicePerSector,fChebZSlicePerSide,1.0f);
3823 0 : fChebCorr->SetUseFloatPrec(kFALSE);
3824 0 : fChebCorr->SetRun(fRun);
3825 0 : fChebCorr->SetTimeStampStart(fTMin);
3826 0 : fChebCorr->SetTimeStampEnd(fTMax);
3827 0 : fChebCorr->SetTimeDependent(kFALSE);
3828 0 : fChebCorr->SetUseZ2R(kTRUE);
3829 : //
3830 0 : if (fBz> 0.01) fChebCorr->SetFieldType(AliTPCChebCorr::kFieldPos);
3831 0 : else if (fBz<-0.01) fChebCorr->SetFieldType(AliTPCChebCorr::kFieldNeg);
3832 0 : else fChebCorr->SetFieldType(AliTPCChebCorr::kFieldZero);
3833 : // Note: to create universal map, set manually SetFieldType(AliTPCChebCorr::kFieldAny)
3834 :
3835 0 : SetUsedInstance(this);
3836 0 : int npCheb[kResDim][2];
3837 0 : for (int i=0;i<kResDim;i++) { // do we need to determine N nodes automatically?
3838 0 : int nbFauto = TMath::Max(int(fNY2XBins*1.2),fNY2XBins+3);
3839 0 : int nbZauto = TMath::Max(int(fNZ2XBins*1.2),fNZ2XBins+3);
3840 0 : npCheb[i][0] = fNPCheb[i][0];
3841 0 : npCheb[i][1] = fNPCheb[i][1];
3842 0 : if (npCheb[i][0]<1) {
3843 0 : npCheb[i][0] = nbFauto; // 1st dimension: sector coordinate y/x
3844 0 : AliInfoF("Nnodes for Cheb.%4s segmentation in %4s is set to %2d",kResName[i],kVoxName[kVoxF],nbFauto);
3845 0 : }
3846 0 : if (npCheb[i][1]<1) {
3847 0 : npCheb[i][1] = nbZauto; // 2nd dimension: z/x
3848 0 : AliInfoF("Nnodes for Cheb.%4s segmentation in %4s is set to %2d",kResName[i],kVoxName[kVoxZ],nbZauto);
3849 0 : }
3850 : }
3851 0 : fChebCorr->Parameterize(trainCorr,kResDim,npCheb,fChebPrecD);
3852 : //
3853 : // register tracks rate for lumi weighting
3854 0 : fChebCorr->SetTracksRate(ExtractTrackRate());
3855 : //
3856 : // calculate weighted lumi
3857 0 : fLumiCTPGraph = AliLumiTools::GetLumiFromCTP(fRun);
3858 0 : fLumiCOG = fChebCorr->GetLuminosityCOG(fLumiCTPGraph);
3859 : //
3860 0 : AliSysInfo::AddStamp("CreateCorrectionObject",1,0,0,0);
3861 0 : }
3862 :
3863 : //________________________________________________________________
3864 : TH1F* AliTPCDcalibRes::ExtractTrackRate() const
3865 : {
3866 : // create histo with used tracks per timestamp
3867 : TH1F* hTr = 0;
3868 0 : if (fTracksRate) {
3869 0 : const float *tarr = fTracksRate->GetArray();
3870 0 : int nb = fTracksRate->GetNbinsX();
3871 : int bin0=1,bin1=nb;
3872 0 : while(tarr[bin0]<0.5 && bin0<=nb) bin0++; // find first significant time bin
3873 0 : while(tarr[bin1]<0.5 && bin1>=bin0) bin1--; // find last significant time bin
3874 0 : nb = bin1 - bin0 + 1;
3875 0 : Long64_t tmn = (Long64_t)fTracksRate->GetBinCenter(bin0);
3876 0 : Long64_t tmx = (Long64_t)fTracksRate->GetBinCenter(bin1);
3877 0 : hTr = new TH1F(Form("TrackRate%lld_%lld",tmn,tmx),"TracksRate",nb,
3878 0 : fTracksRate->GetBinLowEdge(bin0),fTracksRate->GetBinLowEdge(bin1+1));
3879 0 : hTr->SetDirectory(0);
3880 0 : for (int ib=0;ib<nb;ib++) {
3881 0 : int ibh = ib+bin0;
3882 0 : hTr->SetBinContent(ib+1,tarr[ibh]);
3883 : }
3884 0 : }
3885 0 : else AliError("TracksRate accumulation histo was not initialized");
3886 0 : return hTr;
3887 0 : }
3888 :
3889 : //________________________________________________________________
3890 : void AliTPCDcalibRes::InitBinning()
3891 : {
3892 : // initialize binning structures
3893 : //
3894 : // X binning
3895 0 : if (fNXBins>0 && fNXBins<kNPadRows) {
3896 0 : AliInfoF("X-binning: uniform %d bins from %.2f to %.2f",fNXBins,kMinX,kMaxX);
3897 0 : fDXI = fNXBins/(kMaxX-kMinX);
3898 0 : fDX = 1.0f/fDXI;
3899 0 : fUniformBins[kVoxX] = kTRUE;
3900 0 : }
3901 : else {
3902 0 : fNXBins = kNPadRows;
3903 0 : AliInfo("X-binning: bin per pad-row");
3904 0 : fUniformBins[kVoxX] = kFALSE;
3905 0 : fDX = kTPCRowDX[0];
3906 0 : fDXI = 1.f/fDX; // should not be used
3907 : }
3908 : //
3909 : // Y binning
3910 0 : if (fNXBins<1) AliFatal("X bins must be initialized first");
3911 0 : fMaxY2X = new Float_t[fNXBins]; // max Y/X at each X bin, account for dead zones
3912 0 : fDY2XI = new Float_t[fNXBins]; // inverse of Y/X bin size at given X bin
3913 0 : fDY2X = new Float_t[fNXBins]; // Y/X bin size at given X bin
3914 : //
3915 0 : const float kMaxY2X = TMath::Tan(0.5f*kSecDPhi);
3916 :
3917 0 : for (int ix=0;ix<fNXBins;ix++) {
3918 0 : float x = GetX(ix);
3919 0 : fMaxY2X[ix] = kMaxY2X - kDeadZone/x;
3920 0 : fDY2XI[ix] = fNY2XBins / (2.f*fMaxY2X[ix]);
3921 0 : fDY2X[ix] = 1.f/fDY2XI[ix];
3922 : }
3923 : //
3924 0 : fDZ2XI = fNZ2XBins/kMaxZ2X;
3925 0 : fDZ2X = 1.0f/fDZ2XI;
3926 : //
3927 0 : fNBins[kVoxX] = fNXBins;
3928 0 : fNBins[kVoxF] = fNY2XBins;
3929 0 : fNBins[kVoxZ] = fNZ2XBins;
3930 :
3931 0 : fNGVoxPerSector = fNY2XBins*fNZ2XBins*fNXBins;
3932 :
3933 0 : }
3934 :
3935 : //________________________________________________________________
3936 : Int_t AliTPCDcalibRes::GetXBin(float x)
3937 : {
3938 : // convert X to bin ID, following pad row widths
3939 0 : if (fUniformBins[kVoxX]) {
3940 0 : int ix = (x-kMinX)*fDXI;
3941 0 : if (ix<0) return 0;
3942 0 : else if (ix>=fNXBins) return fNXBins-1;
3943 0 : return ix;
3944 : }
3945 : else {
3946 : int ix;
3947 0 : if (x<kTPCRowX[kNRowIROC-1]+0.5*kTPCRowDX[kNRowIROC-1]) { // uniform pad size in IROC
3948 0 : ix = (x-(kTPCRowX[0]-kTPCRowDX[0]*0.5))/kTPCRowDX[0];
3949 0 : if (ix<0) ix = 0;
3950 0 : }
3951 : // uniform pad size in OROC2
3952 0 : else if ( x>= kTPCRowX[kNRowIROC+kNRowOROC1]-0.5*kTPCRowDX[kNRowIROC+kNRowOROC1] ) {
3953 0 : ix = (x-(kTPCRowX[kNRowIROC+kNRowOROC1]-0.5*kTPCRowDX[kNRowIROC+kNRowOROC1]))/kTPCRowDX[kNPadRows-1]
3954 0 : + kNRowIROC + kNRowOROC1;
3955 0 : if (ix>=kNPadRows) ix = kNPadRows-1;
3956 0 : }
3957 : else { // uniform pad size in OROC1
3958 0 : ix = (x-(kTPCRowX[kNRowIROC]-0.5*kTPCRowDX[kNRowIROC]))/kTPCRowDX[kNRowIROC] + kNRowIROC;
3959 0 : if (ix<kNRowIROC) { // we go between IROC and OROC?
3960 0 : if (x> 0.5*(kTPCRowX[kNRowIROC-1]+kTPCRowX[kNRowIROC])) ix = kNRowIROC; // 1st OROC1 row
3961 : else ix = kNRowIROC-1;
3962 : }
3963 :
3964 : }
3965 : return ix;
3966 : //
3967 : }
3968 0 : }
3969 :
3970 : //________________________________________________________________
3971 : Int_t AliTPCDcalibRes::GetRowID(float x)
3972 : {
3973 : // return row ID
3974 : int ix;
3975 0 : if (x<kTPCRowX[kNRowIROC-1]+0.5*kTPCRowDX[kNRowIROC-1]) { // uniform pad size in IROC
3976 0 : ix = (x-(kTPCRowX[0]-kTPCRowDX[0]*0.5))/kTPCRowDX[0];
3977 0 : if (ix<0) ix = -1;
3978 0 : }
3979 : // uniform pad size in OROC2
3980 0 : else if ( x>= kTPCRowX[kNRowIROC+kNRowOROC1]-0.5*kTPCRowDX[kNRowIROC+kNRowOROC1] ) {
3981 0 : ix = (x-(kTPCRowX[kNRowIROC+kNRowOROC1]-0.5*kTPCRowDX[kNRowIROC+kNRowOROC1]))/kTPCRowDX[kNPadRows-1]
3982 0 : + kNRowIROC + kNRowOROC1;
3983 0 : if (ix>=kNPadRows) ix = -2;
3984 0 : }
3985 : else { // uniform pad size in OROC1
3986 0 : ix = (x-(kTPCRowX[kNRowIROC]-0.5*kTPCRowDX[kNRowIROC]))/kTPCRowDX[kNRowIROC] + kNRowIROC;
3987 0 : if (ix<kNRowIROC) { // we go between IROC and OROC?
3988 : ix = -3;
3989 : }
3990 :
3991 : }
3992 0 : return ix;
3993 : }
3994 :
3995 : //_____________________________________________________
3996 : Bool_t AliTPCDcalibRes::FindVoxelBin(int sectID, float x, float y, float z, UChar_t bin[kVoxHDim],float voxVars[kVoxHDim])
3997 : {
3998 : // define voxel variables and bin
3999 : //
4000 : //
4001 : // track Z/X bin
4002 0 : voxVars[kVoxZ] = z/x;
4003 0 : if (TMath::Abs(voxVars[kVoxZ])>kMaxZ2X) return kFALSE;
4004 0 : int binZ = GetZ2XBinExact( sectID<kNSect ? voxVars[kVoxZ] : -voxVars[kVoxZ] );
4005 0 : if (binZ<0) return kFALSE;
4006 0 : bin[kVoxZ] = binZ;
4007 : //
4008 : // track X bin
4009 0 : voxVars[kVoxX] = x;
4010 0 : int binX = GetXBinExact(x); // exact matching to pad-rows
4011 0 : if (binX<0 || binX>=fNXBins) return kFALSE;
4012 0 : bin[kVoxX] = binX;
4013 : //
4014 : // track Y/X bin accounting for sector edges dead zones
4015 0 : voxVars[kVoxF] = y/x;
4016 0 : int binF = GetY2XBinExact(voxVars[kVoxF],binX);
4017 0 : if (binF<0||binF>=fNY2XBins) return kFALSE;
4018 0 : bin[kVoxF] = binF;
4019 : //
4020 0 : return kTRUE;
4021 0 : }
4022 :
4023 : //=====================================================================
4024 : ///////////////////////////////////////////////////
4025 : //
4026 : // This is a special non-meber f-n for cheb. learning
4027 : //
4028 : void trainCorr(int row, float* tzLoc, float* corrLoc)
4029 : {
4030 : // Cheb. object training f-n: compute correction for the point
4031 : //
4032 : // xtzLoc: y2x and z2x in sector frame
4033 : // corrLoc: x, y, z corrections in sector frame
4034 : // when called with pointer at 0, row will set the sector/side
4035 : // (row should be sector in 0-35 or 0-71 format)
4036 : static int sector=0;
4037 0 : if (!tzLoc || !corrLoc) {
4038 0 : sector = row%AliTPCDcalibRes::kNSect2;
4039 0 : printf("training Sector%d\n",sector);
4040 0 : return;
4041 : }
4042 : //
4043 0 : AliTPCDcalibRes* calib = AliTPCDcalibRes::GetUsedInstance();
4044 :
4045 0 : float x = AliTPCDcalibRes::GetTPCRowX(row);
4046 0 : float dist[AliTPCDcalibRes::kResDim] = {0};
4047 0 : int xbin = calib->GetNXBins()==AliTPCDcalibRes::kNPadRows ? row : calib->GetXBin(x);
4048 0 : if (calib->GetXBinIgnored(sector,xbin)) return;
4049 0 : float y2x = tzLoc[0];
4050 0 : float z2x = tzLoc[1];
4051 : //
4052 0 : Bool_t res = calib->GetSmoothEstimate(sector, x, y2x, z2x, 0xff, dist);
4053 0 : if (!res) { printf("Failed to evaluate smooth distortion\n"); exit(1); }
4054 :
4055 : /*
4056 : // Marian stored Z track coordinate instead of cluster one, need to correct for this
4057 : if (fApplyZt2Zc) {
4058 : float deriv[AliTPCDcalibRes::kResDim*3];
4059 : const double inversionEps = 20e-4; // when inverting, stop Newton-Raphson iterations at this eps
4060 : const int inversionMaxIt = 3; // when inverting, stop Newton-Raphson after some numbers of iterations
4061 : double change = 0, xInv = 1./x;
4062 : float zc = z2x*x;
4063 : float zt = zc + dist[kResZ]; // 1st guess on true zt at measured zc
4064 : //
4065 : // use Newton-Raphson method for NDLocal inversion to get zt = F(zc) from
4066 : // dz == zt - zc = NDLoc(zt) -> zc - [zt - NDLoc(zt)]=0
4067 : // -> zt_{i+1} = zt_i - (zc - [zt - NDLoct(zt_i)])/(-d[zt_i-NDLoc(zt_i)]/dz)
4068 : //
4069 : int it = 0;
4070 : do {
4071 : z2x = zt*xInv;
4072 : res = GetSmoothEstimate(sector, x, y2x, z2x, dist, deriv);
4073 : if (!res) {printf("Failed to evaluate smooth distortion\n");exit(1);}
4074 : double bot = 1. - deriv[kResZ*3+2]*xInv; // dF(zt_i)/dz
4075 : if (TMath::Abs(bot)<1e-6) break;
4076 : double top = zc - (zt - dist[kResZ]); // zc - F(zt_i)
4077 : double change = top/bot;
4078 : // printf("It %d Eps:%+e, Zc:%+e, Zt:%+e Ztn:%+e DZ:%+e | dztmp: :%+e DD:%+e\n",
4079 : // it,change,zc,zt,zt+change,dz,dztmp,deriv[kndZ2R]*xInv);
4080 : zt += change;
4081 : it++;
4082 : } while(it<inversionMaxIt && TMath::Abs(change)>inversionEps);
4083 : //
4084 : // now query at fixed Z2X
4085 : z2x = zt*xInv;
4086 : res = GetSmoothEstimate(sector, x, y2x, z2x, dist);
4087 : if (!res) {printf("Failed to evaluate smooth distortion\n");exit(1);}
4088 : }
4089 : */
4090 :
4091 0 : for (int i=0;i<AliTPCDcalibRes::kResDim;i++) corrLoc[i] = dist[i];
4092 : //
4093 0 : }
4094 : //======================================================================================
|
