Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 :
17 : /// \class AliTPCcalibDButil
18 : /// \brief Class providing the calculation of derived quantities (mean,rms,fits,...) of calibration entries
19 :
20 : #include <TMath.h>
21 : #include <TVectorT.h>
22 : #include <TObjArray.h>
23 : #include <TGraph.h>
24 : #include <TFile.h>
25 : #include <TDirectory.h>
26 : #include <TMap.h>
27 : #include <TGraphErrors.h>
28 : #include <AliCDBStorage.h>
29 : #include <AliDCSSensorArray.h>
30 : #include <AliTPCSensorTempArray.h>
31 : #include <AliDCSSensor.h>
32 : #include <AliLog.h>
33 : #include <AliCDBEntry.h>
34 : #include <AliCDBManager.h>
35 : #include <AliCDBId.h>
36 : #include <AliSplineFit.h>
37 : #include "AliTPCcalibDB.h"
38 : #include "AliTPCCalPad.h"
39 : #include "AliTPCCalROC.h"
40 : #include "AliTPCROC.h"
41 : #include "AliTPCmapper.h"
42 : #include "AliTPCParam.h"
43 : #include "AliTPCCalibRaw.h"
44 : #include "AliTPCPreprocessorOnline.h"
45 : #include "AliTPCdataQA.h"
46 : #include "AliLog.h"
47 : #include "AliTPCcalibDButil.h"
48 : #include "AliTPCCalibVdrift.h"
49 : #include "AliMathBase.h"
50 : #include "AliRelAlignerKalman.h"
51 : #include "TTree.h"
52 : #include "TROOT.h"
53 : #include "TKey.h"
54 :
55 : const Float_t kAlmost0=1.e-30;
56 :
57 : /// \cond CLASSIMP
58 24 : ClassImp(AliTPCcalibDButil)
59 : /// \endcond
60 : AliTPCcalibDButil::AliTPCcalibDButil() :
61 3 : TObject(),
62 3 : fCalibDB(0),
63 3 : fPadNoise(0x0),
64 3 : fPedestals(0x0),
65 3 : fPulserTmean(0x0),
66 3 : fPulserTrms(0x0),
67 3 : fPulserQmean(0x0),
68 9 : fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
69 3 : fCETmean(0x0),
70 3 : fCETrms(0x0),
71 3 : fCEQmean(0x0),
72 3 : fALTROMasked(0x0),
73 3 : fCalibRaw(0x0),
74 3 : fDataQA(0x0),
75 3 : fRefMap(0x0),
76 3 : fCurrentRefMap(0x0),
77 3 : fRefValidity(""),
78 3 : fRefPadNoise(0x0),
79 3 : fRefPedestals(0x0),
80 3 : fRefPedestalMasked(0x0),
81 3 : fRefPulserTmean(0x0),
82 3 : fRefPulserTrms(0x0),
83 3 : fRefPulserQmean(0x0),
84 9 : fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
85 3 : fRefPulserMasked(0x0),
86 3 : fRefCETmean(0x0),
87 3 : fRefCETrms(0x0),
88 3 : fRefCEQmean(0x0),
89 3 : fRefCEMasked(0x0),
90 3 : fRefALTROFPED(0x0),
91 3 : fRefALTROZsThr(0x0),
92 3 : fRefALTROAcqStart(0x0),
93 3 : fRefALTROAcqStop(0x0),
94 3 : fRefALTROMasked(0x0),
95 3 : fRefCalibRaw(0x0),
96 3 : fRefDataQA(0x0),
97 3 : fGoofieArray(0x0),
98 9 : fMapper(new AliTPCmapper(0x0)),
99 3 : fNpulserOutliers(-1),
100 3 : fIrocTimeOffset(0),
101 3 : fCETmaxLimitAbs(1.5),
102 3 : fPulTmaxLimitAbs(1.5),
103 3 : fPulQmaxLimitAbs(5),
104 3 : fPulQminLimit(11),
105 3 : fRuns(0), // run list with OCDB info
106 3 : fRunsStart(0), // start time for given run
107 3 : fRunsStop(0) // stop time for given run
108 15 : {
109 : //
110 : // Default ctor
111 : //
112 6 : }
113 : //_____________________________________________________________________________________
114 : AliTPCcalibDButil::~AliTPCcalibDButil()
115 0 : {
116 : /// dtor
117 :
118 0 : delete fPulserOutlier;
119 0 : delete fRefPulserOutlier;
120 0 : delete fMapper;
121 0 : if (fRefPadNoise) delete fRefPadNoise;
122 0 : if (fRefPedestals) delete fRefPedestals;
123 0 : if (fRefPedestalMasked) delete fRefPedestalMasked;
124 0 : if (fRefPulserTmean) delete fRefPulserTmean;
125 0 : if (fRefPulserTrms) delete fRefPulserTrms;
126 0 : if (fRefPulserQmean) delete fRefPulserQmean;
127 0 : if (fRefPulserMasked) delete fRefPulserMasked;
128 0 : if (fRefCETmean) delete fRefCETmean;
129 0 : if (fRefCETrms) delete fRefCETrms;
130 0 : if (fRefCEQmean) delete fRefCEQmean;
131 0 : if (fRefCEMasked) delete fRefCEMasked;
132 0 : if (fRefALTROFPED) delete fRefALTROFPED;
133 0 : if (fRefALTROZsThr) delete fRefALTROZsThr;
134 0 : if (fRefALTROAcqStart) delete fRefALTROAcqStart;
135 0 : if (fRefALTROAcqStop) delete fRefALTROAcqStop;
136 0 : if (fRefALTROMasked) delete fRefALTROMasked;
137 0 : if (fRefCalibRaw) delete fRefCalibRaw;
138 0 : if (fCurrentRefMap) delete fCurrentRefMap;
139 0 : }
140 : //_____________________________________________________________________________________
141 : void AliTPCcalibDButil::UpdateFromCalibDB()
142 : {
143 : /// Update pointers from calibDB
144 :
145 0 : if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
146 0 : fCalibDB->UpdateNonRec(); // load all infromation now
147 0 : fPadNoise=fCalibDB->GetPadNoise();
148 0 : fPedestals=fCalibDB->GetPedestals();
149 0 : fPulserTmean=fCalibDB->GetPulserTmean();
150 0 : fPulserTrms=fCalibDB->GetPulserTrms();
151 0 : fPulserQmean=fCalibDB->GetPulserQmean();
152 0 : fCETmean=fCalibDB->GetCETmean();
153 0 : fCETrms=fCalibDB->GetCETrms();
154 0 : fCEQmean=fCalibDB->GetCEQmean();
155 0 : fALTROMasked=fCalibDB->GetALTROMasked();
156 0 : fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
157 0 : fCalibRaw=fCalibDB->GetCalibRaw();
158 0 : fDataQA=fCalibDB->GetDataQA();
159 0 : UpdatePulserOutlierMap();
160 : // SetReferenceRun();
161 : // UpdateRefDataFromOCDB();
162 0 : }
163 : //_____________________________________________________________________________________
164 : void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
165 : Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
166 : {
167 : /// Process the CE data for this run
168 : /// the return TVectorD arrays contian the results of the fit
169 : /// noutliersCE contains the number of pads marked as outliers,
170 : /// not including masked and edge pads
171 :
172 : //retrieve CE and ALTRO data
173 0 : if (!fCETmean){
174 0 : TString fitString(fitFormula);
175 0 : fitString.ReplaceAll("++","#");
176 0 : Int_t ndim=fitString.CountChar('#')+2;
177 0 : fitResultsA.ResizeTo(ndim);
178 0 : fitResultsC.ResizeTo(ndim);
179 0 : fitResultsA.Zero();
180 0 : fitResultsC.Zero();
181 0 : noutliersCE=-1;
182 : return;
183 0 : }
184 0 : noutliersCE=0;
185 : //create outlier map
186 : AliTPCCalPad *out=0;
187 0 : if (outCE) out=outCE;
188 0 : else out=new AliTPCCalPad("outCE","outCE");
189 : AliTPCCalROC *rocMasked=0x0;
190 : //loop over all channels
191 0 : for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
192 0 : AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
193 0 : if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
194 0 : AliTPCCalROC *rocOut=out->GetCalROC(iroc);
195 0 : if (!rocData) {
196 0 : noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
197 0 : rocOut->Add(1.);
198 0 : continue;
199 : }
200 : //add time offset to IROCs
201 0 : if (iroc<AliTPCROC::Instance()->GetNInnerSector())
202 0 : rocData->Add(fIrocTimeOffset);
203 : //select outliers
204 0 : UInt_t nrows=rocData->GetNrows();
205 0 : for (UInt_t irow=0;irow<nrows;++irow){
206 0 : UInt_t npads=rocData->GetNPads(irow);
207 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
208 0 : rocOut->SetValue(irow,ipad,0);
209 : //exclude masked pads
210 0 : if (rocMasked && rocMasked->GetValue(irow,ipad)) {
211 0 : rocOut->SetValue(irow,ipad,1);
212 0 : continue;
213 : }
214 : //exclude first two rows in IROC and last two rows in OROC
215 0 : if (iroc<36){
216 0 : if (irow<2) rocOut->SetValue(irow,ipad,1);
217 : } else {
218 0 : if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
219 : }
220 : //exclude edge pads
221 0 : if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
222 0 : Float_t valTmean=rocData->GetValue(irow,ipad);
223 : //exclude values that are exactly 0
224 0 : if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
225 0 : rocOut->SetValue(irow,ipad,1);
226 0 : ++noutliersCE;
227 0 : }
228 : // exclude channels with too large variations
229 0 : if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
230 0 : rocOut->SetValue(irow,ipad,1);
231 0 : ++noutliersCE;
232 0 : }
233 0 : }
234 : }
235 0 : }
236 : //perform fit
237 0 : TMatrixD dummy;
238 0 : Float_t chi2Af,chi2Cf;
239 0 : fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
240 0 : chi2A=chi2Af;
241 0 : chi2C=chi2Cf;
242 0 : if (!outCE) delete out;
243 0 : }
244 : //_____________________________________________________________________________________
245 : void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
246 : TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
247 : Float_t &driftTimeA, Float_t &driftTimeC )
248 : {
249 : /// Calculate statistical information from the CE graphs for drift time and charge
250 :
251 : //reset arrays
252 0 : vecTEntries.ResizeTo(72);
253 0 : vecTMean.ResizeTo(72);
254 0 : vecTRMS.ResizeTo(72);
255 0 : vecTMedian.ResizeTo(72);
256 0 : vecQEntries.ResizeTo(72);
257 0 : vecQMean.ResizeTo(72);
258 0 : vecQRMS.ResizeTo(72);
259 0 : vecQMedian.ResizeTo(72);
260 0 : vecTEntries.Zero();
261 0 : vecTMean.Zero();
262 0 : vecTRMS.Zero();
263 0 : vecTMedian.Zero();
264 0 : vecQEntries.Zero();
265 0 : vecQMean.Zero();
266 0 : vecQRMS.Zero();
267 0 : vecQMedian.Zero();
268 0 : driftTimeA=0;
269 0 : driftTimeC=0;
270 0 : TObjArray *arrT=fCalibDB->GetCErocTtime();
271 0 : TObjArray *arrQ=fCalibDB->GetCErocQtime();
272 0 : if (arrT){
273 0 : for (Int_t isec=0;isec<74;++isec){
274 0 : TGraph *gr=(TGraph*)arrT->At(isec);
275 0 : if (!gr) continue;
276 0 : TVectorD values;
277 0 : Int_t npoints = gr->GetN();
278 0 : values.ResizeTo(npoints);
279 : Int_t nused =0;
280 : //skip first points, theres always some problems with finding the CE position
281 0 : for (Int_t ipoint=4; ipoint<npoints; ipoint++){
282 0 : if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
283 0 : values[nused]=gr->GetY()[ipoint];
284 0 : nused++;
285 0 : }
286 : }
287 : //
288 0 : if (isec<72) vecTEntries[isec]= nused;
289 0 : if (nused>1){
290 0 : if (isec<72){
291 0 : vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
292 0 : vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
293 0 : vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
294 0 : } else if (isec==72){
295 0 : driftTimeA=TMath::Median(nused,values.GetMatrixArray());
296 0 : } else if (isec==73){
297 0 : driftTimeC=TMath::Median(nused,values.GetMatrixArray());
298 0 : }
299 : }
300 0 : }
301 0 : }
302 0 : if (arrQ){
303 0 : for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
304 0 : TGraph *gr=(TGraph*)arrQ->At(isec);
305 0 : if (!gr) continue;
306 0 : TVectorD values;
307 0 : Int_t npoints = gr->GetN();
308 0 : values.ResizeTo(npoints);
309 : Int_t nused =0;
310 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
311 0 : if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
312 0 : values[nused]=gr->GetY()[ipoint];
313 0 : nused++;
314 0 : }
315 : }
316 : //
317 0 : vecQEntries[isec]= nused;
318 0 : if (nused>1){
319 0 : vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
320 0 : vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
321 0 : vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
322 0 : }
323 0 : }
324 0 : }
325 0 : }
326 :
327 : //_____________________________________________________________________________________
328 : void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
329 : TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
330 : Int_t &nonMaskedZero, Int_t &nNaN)
331 : {
332 : /// process noise data
333 : /// vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
334 : /// OROCs small pads [2] and OROCs large pads [3]
335 : /// vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
336 : /// nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
337 :
338 : //set proper size and reset
339 : const UInt_t infoSize=4;
340 0 : vNoiseMean.ResizeTo(infoSize);
341 0 : vNoiseMeanSenRegions.ResizeTo(infoSize);
342 0 : vNoiseRMS.ResizeTo(infoSize);
343 0 : vNoiseRMSSenRegions.ResizeTo(infoSize);
344 0 : vNoiseMean.Zero();
345 0 : vNoiseMeanSenRegions.Zero();
346 0 : vNoiseRMS.Zero();
347 0 : vNoiseRMSSenRegions.Zero();
348 0 : nonMaskedZero=0;
349 0 : nNaN=0;
350 : //counters
351 0 : TVectorD c(infoSize);
352 0 : TVectorD cs(infoSize);
353 : //tpc parameters
354 0 : AliTPCParam par;
355 0 : par.Update();
356 : //retrieve noise and ALTRO data
357 0 : if (!fPadNoise) return;
358 : AliTPCCalROC *rocMasked=0x0;
359 : //create IROC, OROC1, OROC2 and sensitive region masks
360 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
361 0 : AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
362 0 : if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
363 0 : UInt_t nrows=noiseROC->GetNrows();
364 0 : for (UInt_t irow=0;irow<nrows;++irow){
365 0 : UInt_t npads=noiseROC->GetNPads(irow);
366 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
367 : //don't use masked channels;
368 0 : if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
369 0 : Float_t noiseVal=noiseROC->GetValue(irow,ipad);
370 : //check if noise==0
371 0 : if (noiseVal<kAlmost0) {
372 0 : ++nonMaskedZero;
373 0 : continue;
374 : }
375 : //check for nan
376 0 : if ( !(noiseVal<10000000) ){
377 0 : AliInfo(Form("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal));
378 0 : ++nNaN;
379 0 : continue;
380 : }
381 0 : Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
382 : Int_t masksen=1; // sensitive pards are not masked (0)
383 0 : if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
384 0 : if (isec<AliTPCROC::Instance()->GetNInnerSector()){
385 : //IROCs
386 0 : if (irow>19&&irow<46){
387 0 : if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
388 : }
389 : Int_t type=1;
390 0 : vNoiseMean[type]+=noiseVal;
391 0 : vNoiseRMS[type]+=noiseVal*noiseVal;
392 0 : ++c[type];
393 0 : if (!masksen){
394 0 : vNoiseMeanSenRegions[type]+=noiseVal;
395 0 : vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
396 0 : ++cs[type];
397 0 : }
398 0 : } else {
399 : //OROCs
400 : //define sensive regions
401 0 : if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
402 0 : if ( irow>75 ){
403 0 : Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
404 0 : if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
405 0 : }
406 0 : if ((Int_t)irow<par.GetNRowUp1()){
407 : //OROC1
408 : Int_t type=2;
409 0 : vNoiseMean[type]+=noiseVal;
410 0 : vNoiseRMS[type]+=noiseVal*noiseVal;
411 0 : ++c[type];
412 0 : if (!masksen){
413 0 : vNoiseMeanSenRegions[type]+=noiseVal;
414 0 : vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
415 0 : ++cs[type];
416 0 : }
417 0 : }else{
418 : //OROC2
419 : Int_t type=3;
420 0 : vNoiseMean[type]+=noiseVal;
421 0 : vNoiseRMS[type]+=noiseVal*noiseVal;
422 0 : ++c[type];
423 0 : if (!masksen){
424 0 : vNoiseMeanSenRegions[type]+=noiseVal;
425 0 : vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
426 0 : ++cs[type];
427 0 : }
428 : }
429 : }
430 : //whole tpc
431 : Int_t type=0;
432 0 : vNoiseMean[type]+=noiseVal;
433 0 : vNoiseRMS[type]+=noiseVal*noiseVal;
434 0 : ++c[type];
435 0 : if (!masksen){
436 0 : vNoiseMeanSenRegions[type]+=noiseVal;
437 0 : vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
438 0 : ++cs[type];
439 0 : }
440 0 : }//end loop pads
441 : }//end loop rows
442 : }//end loop sectors (rocs)
443 :
444 : //calculate mean and RMS
445 : const Double_t verySmall=0.0000000001;
446 0 : for (UInt_t i=0;i<infoSize;++i){
447 : Double_t mean=0;
448 : Double_t rms=0;
449 : Double_t meanSen=0;
450 : Double_t rmsSen=0;
451 :
452 0 : if (c[i]>verySmall){
453 0 : AliInfo(Form("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]));
454 0 : mean=vNoiseMean[i]/c[i];
455 0 : rms=vNoiseRMS[i];
456 0 : rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
457 0 : }
458 0 : vNoiseMean[i]=mean;
459 0 : vNoiseRMS[i]=rms;
460 :
461 0 : if (cs[i]>verySmall){
462 0 : meanSen=vNoiseMeanSenRegions[i]/cs[i];
463 0 : rmsSen=vNoiseRMSSenRegions[i];
464 0 : rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
465 0 : }
466 0 : vNoiseMeanSenRegions[i]=meanSen;
467 0 : vNoiseRMSSenRegions[i]=rmsSen;
468 : }
469 0 : }
470 :
471 : //_____________________________________________________________________________________
472 : void AliTPCcalibDButil::ProcessQAData(TVectorD &vQaOcc, TVectorD &vQaQtot,
473 : TVectorD &vQaQmax)
474 : {
475 : /// process QA data
476 : ///
477 : /// vQaOcc/Qtot/Qmax contains the Mean occupancy/Qtot/Qmax for each sector
478 :
479 :
480 : const UInt_t infoSize = 72;
481 : //reset counters to error number
482 0 : vQaOcc.ResizeTo(infoSize);
483 0 : vQaOcc.Zero();
484 0 : vQaQtot.ResizeTo(infoSize);
485 0 : vQaQtot.Zero();
486 0 : vQaQmax.ResizeTo(infoSize);
487 0 : vQaQmax.Zero();
488 : //counter
489 : //retrieve pulser and ALTRO data
490 :
491 0 : if (!fDataQA) {
492 :
493 0 : AliInfo("No QA data");
494 0 : return;
495 : }
496 0 : if (fDataQA->GetEventCounter()<=0) {
497 :
498 0 : AliInfo("No QA data");
499 0 : return; // no data processed
500 : }
501 : //
502 0 : fDataQA->Analyse();
503 :
504 0 : TVectorD normOcc(infoSize);
505 0 : TVectorD normQ(infoSize);
506 :
507 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
508 :
509 0 : AliInfo(Form("Sector %d\n", isec));
510 0 : AliTPCCalROC* occupancyROC = fDataQA->GetNoThreshold()->GetCalROC(isec);
511 0 : AliTPCCalROC* nclusterROC = fDataQA->GetNLocalMaxima()->GetCalROC(isec);
512 0 : AliTPCCalROC* qROC = fDataQA->GetMeanCharge()->GetCalROC(isec);
513 0 : AliTPCCalROC* qmaxROC = fDataQA->GetMaxCharge()->GetCalROC(isec);
514 0 : if (!occupancyROC) continue;
515 0 : if (!nclusterROC) continue;
516 0 : if (!qROC) continue;
517 0 : if (!qmaxROC) continue;
518 :
519 0 : const UInt_t nchannels=occupancyROC->GetNchannels();
520 :
521 0 : AliInfo(Form("Nchannels %d\n", nchannels));
522 :
523 0 : for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
524 :
525 0 : vQaOcc[isec] += occupancyROC->GetValue(ichannel);
526 0 : ++normOcc[isec];
527 :
528 0 : Float_t nClusters = nclusterROC->GetValue(ichannel);
529 0 : normQ[isec] += nClusters;
530 0 : vQaQtot[isec]+=nClusters*qROC->GetValue(ichannel);
531 0 : vQaQmax[isec]+=nClusters*qmaxROC->GetValue(ichannel);
532 : }
533 0 : }
534 :
535 : //calculate mean values
536 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
537 :
538 0 : if (normOcc[isec]>0) vQaOcc[isec] /= normOcc[isec];
539 0 : else vQaOcc[isec] = 0;
540 :
541 0 : if (normQ[isec]>0) {
542 0 : vQaQtot[isec] /= normQ[isec];
543 0 : vQaQmax[isec] /= normQ[isec];
544 0 : }else {
545 :
546 0 : vQaQtot[isec] = 0;
547 0 : vQaQmax[isec] = 0;
548 : }
549 : }
550 0 : }
551 :
552 : //_____________________________________________________________________________________
553 : void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
554 : {
555 : /// Process the Pulser information
556 : /// vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
557 :
558 : const UInt_t infoSize=4;
559 : //reset counters to error number
560 0 : vMeanTime.ResizeTo(infoSize);
561 0 : vMeanTime.Zero();
562 : //counter
563 0 : TVectorD c(infoSize);
564 : //retrieve pulser and ALTRO data
565 0 : if (!fPulserTmean) return;
566 : //
567 : //get Outliers
568 : AliTPCCalROC *rocOut=0x0;
569 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
570 0 : AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
571 0 : if (!tmeanROC) continue;
572 0 : rocOut=fPulserOutlier->GetCalROC(isec);
573 0 : UInt_t nchannels=tmeanROC->GetNchannels();
574 0 : for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
575 0 : if (rocOut && rocOut->GetValue(ichannel)) continue;
576 0 : Float_t val=tmeanROC->GetValue(ichannel);
577 0 : Int_t type=isec/18;
578 0 : vMeanTime[type]+=val;
579 0 : ++c[type];
580 0 : }
581 0 : }
582 : //calculate mean
583 0 : for (UInt_t itype=0; itype<infoSize; ++itype){
584 0 : if (c[itype]>0) vMeanTime[itype]/=c[itype];
585 0 : else vMeanTime[itype]=0;
586 : }
587 0 : }
588 : //_____________________________________________________________________________________
589 : void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
590 : {
591 : /// Get Values from ALTRO configuration data
592 :
593 0 : nMasked=-1;
594 0 : if (!fALTROMasked) return;
595 0 : nMasked=0;
596 0 : for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
597 0 : AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
598 0 : for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
599 0 : if (rocMasked->GetValue(ichannel)) ++nMasked;
600 : }
601 : }
602 0 : }
603 : //_____________________________________________________________________________________
604 : void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
605 : {
606 : /// Proces Goofie values, return statistical information of the currently set goofieArray
607 : /// The meaning of the entries are given below
608 :
609 : /*
610 : 1 TPC_ANODE_I_A00_STAT
611 : 2 TPC_DVM_CO2
612 : 3 TPC_DVM_DriftVelocity
613 : 4 TPC_DVM_FCageHV
614 : 5 TPC_DVM_GainFar
615 : 6 TPC_DVM_GainNear
616 : 7 TPC_DVM_N2
617 : 8 TPC_DVM_NumberOfSparks
618 : 9 TPC_DVM_PeakAreaFar
619 : 10 TPC_DVM_PeakAreaNear
620 : 11 TPC_DVM_PeakPosFar
621 : 12 TPC_DVM_PeakPosNear
622 : 13 TPC_DVM_PickupHV
623 : 14 TPC_DVM_Pressure
624 : 15 TPC_DVM_T1_Over_P
625 : 16 TPC_DVM_T2_Over_P
626 : 17 TPC_DVM_T_Over_P
627 : 18 TPC_DVM_TemperatureS1
628 : */
629 0 : if (!fGoofieArray){
630 : Int_t nsensors=19;
631 0 : vecEntries.ResizeTo(nsensors);
632 0 : vecMedian.ResizeTo(nsensors);
633 0 : vecMean.ResizeTo(nsensors);
634 0 : vecRMS.ResizeTo(nsensors);
635 0 : vecEntries.Zero();
636 0 : vecMedian.Zero();
637 0 : vecMean.Zero();
638 0 : vecRMS.Zero();
639 : return;
640 : }
641 : Double_t kEpsilon=0.0000000001;
642 : Double_t kBig=100000000000.;
643 0 : Int_t nsensors = fGoofieArray->NumSensors();
644 0 : vecEntries.ResizeTo(nsensors);
645 0 : vecMedian.ResizeTo(nsensors);
646 0 : vecMean.ResizeTo(nsensors);
647 0 : vecRMS.ResizeTo(nsensors);
648 0 : TVectorF values;
649 0 : for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
650 0 : AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
651 0 : if (gsensor && gsensor->GetGraph()){
652 0 : Int_t npoints = gsensor->GetGraph()->GetN();
653 : // filter zeroes
654 0 : values.ResizeTo(npoints);
655 : Int_t nused =0;
656 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
657 0 : if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
658 0 : TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
659 0 : values[nused]=gsensor->GetGraph()->GetY()[ipoint];
660 0 : nused++;
661 0 : }
662 : }
663 : //
664 0 : vecEntries[isensor]= nused;
665 0 : if (nused>1){
666 0 : vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
667 0 : vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
668 0 : vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
669 0 : }
670 0 : }
671 : }
672 0 : }
673 : //_____________________________________________________________________________________
674 : void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
675 : {
676 : /// check the variations of the pedestal data to the reference pedestal data
677 : /// thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
678 :
679 : const Int_t npar=4;
680 0 : TVectorF vThres(npar); //thresholds
681 : Int_t nActive=0; //number of active channels
682 :
683 : //reset and set thresholds
684 0 : pedestalDeviations.ResizeTo(npar);
685 0 : for (Int_t i=0;i<npar;++i){
686 0 : pedestalDeviations.GetMatrixArray()[i]=0;
687 0 : vThres.GetMatrixArray()[i]=(i+1)*.5;
688 : }
689 : //check all needed data is available
690 0 : if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
691 : //loop over all channels
692 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
693 0 : AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
694 0 : AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
695 0 : AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
696 0 : AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
697 0 : UInt_t nrows=mROC->GetNrows();
698 0 : for (UInt_t irow=0;irow<nrows;++irow){
699 0 : UInt_t npads=mROC->GetNPads(irow);
700 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
701 : //don't use masked channels;
702 0 : if (mROC ->GetValue(irow,ipad)) continue;
703 0 : if (mRefROC->GetValue(irow,ipad)) continue;
704 0 : Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
705 0 : for (Int_t i=0;i<npar;++i){
706 0 : if (deviation>vThres[i])
707 0 : ++pedestalDeviations.GetMatrixArray()[i];
708 : }
709 0 : ++nActive;
710 0 : }//end ipad
711 : }//ind irow
712 : }//end isec
713 0 : if (nActive>0){
714 0 : for (Int_t i=0;i<npar;++i){
715 0 : pedestalDeviations.GetMatrixArray()[i]/=nActive;
716 : }
717 0 : }
718 0 : }
719 : //_____________________________________________________________________________________
720 : void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
721 : {
722 : /// check the variations of the noise data to the reference noise data
723 : /// thresholds are 5, 10, 15 and 20 percent respectively.
724 :
725 : const Int_t npar=4;
726 0 : TVectorF vThres(npar); //thresholds
727 : Int_t nActive=0; //number of active channels
728 :
729 : //reset and set thresholds
730 0 : noiseDeviations.ResizeTo(npar);
731 0 : for (Int_t i=0;i<npar;++i){
732 0 : noiseDeviations.GetMatrixArray()[i]=0;
733 0 : vThres.GetMatrixArray()[i]=(i+1)*.05;
734 : }
735 : //check all needed data is available
736 0 : if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
737 : //loop over all channels
738 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
739 0 : AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
740 0 : AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
741 0 : AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
742 0 : AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
743 0 : UInt_t nrows=mROC->GetNrows();
744 0 : for (UInt_t irow=0;irow<nrows;++irow){
745 0 : UInt_t npads=mROC->GetNPads(irow);
746 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
747 : //don't use masked channels;
748 0 : if (mROC ->GetValue(irow,ipad)) continue;
749 0 : if (mRefROC->GetValue(irow,ipad)) continue;
750 0 : if (nRefROC->GetValue(irow,ipad)==0) continue;
751 0 : Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
752 0 : for (Int_t i=0;i<npar;++i){
753 0 : if (deviation>vThres[i])
754 0 : ++noiseDeviations.GetMatrixArray()[i];
755 : }
756 0 : ++nActive;
757 0 : }//end ipad
758 : }//ind irow
759 : }//end isec
760 0 : if (nActive>0){
761 0 : for (Int_t i=0;i<npar;++i){
762 0 : noiseDeviations.GetMatrixArray()[i]/=nActive;
763 : }
764 0 : }
765 0 : }
766 : //_____________________________________________________________________________________
767 : void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
768 : Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
769 : {
770 : /// check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
771 : /// thresholds are .5, 1, 5 and 10 percent respectively.
772 :
773 : const Int_t npar=4;
774 0 : TVectorF vThres(npar); //thresholds
775 : Int_t nActive=0; //number of active channels
776 :
777 : //reset and set thresholds
778 0 : pulserQdeviations.ResizeTo(npar);
779 0 : for (Int_t i=0;i<npar;++i){
780 0 : pulserQdeviations.GetMatrixArray()[i]=0;
781 : }
782 0 : npadsOutOneTB=0;
783 0 : npadsOffAdd=0;
784 0 : varQMean=0;
785 0 : vThres.GetMatrixArray()[0]=.005;
786 0 : vThres.GetMatrixArray()[1]=.01;
787 0 : vThres.GetMatrixArray()[2]=.05;
788 0 : vThres.GetMatrixArray()[3]=.1;
789 : //check all needed data is available
790 0 : if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
791 : //
792 0 : UpdateRefPulserOutlierMap();
793 : //loop over all channels
794 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
795 0 : AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
796 0 : AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
797 0 : AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
798 : // AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
799 0 : AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
800 0 : AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
801 0 : AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
802 0 : Float_t ptmean=ptROC->GetMean(oROC);
803 0 : UInt_t nrows=mROC->GetNrows();
804 0 : for (UInt_t irow=0;irow<nrows;++irow){
805 0 : UInt_t npads=mROC->GetNPads(irow);
806 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
807 : //don't use masked channels;
808 0 : if (mROC ->GetValue(irow,ipad)) continue;
809 0 : if (mRefROC->GetValue(irow,ipad)) continue;
810 : //don't user edge pads
811 0 : if (ipad==0||ipad==npads-1) continue;
812 : //data
813 0 : Float_t pq=pqROC->GetValue(irow,ipad);
814 0 : Float_t pqRef=pqRefROC->GetValue(irow,ipad);
815 0 : Float_t pt=ptROC->GetValue(irow,ipad);
816 : // Float_t ptRef=ptRefROC->GetValue(irow,ipad);
817 : //comparisons q
818 0 : Float_t deviation=TMath::Abs(pqRef)>1e-20?TMath::Abs(pq/pqRef-1):-999;
819 0 : for (Int_t i=0;i<npar;++i){
820 0 : if (deviation>vThres[i])
821 0 : ++pulserQdeviations.GetMatrixArray()[i];
822 : }
823 0 : if (pqRef>11&&pq<11) ++npadsOffAdd;
824 0 : varQMean+=pq-pqRef;
825 : //comparisons t
826 0 : if (TMath::Abs(pt-ptmean)>1) ++npadsOutOneTB;
827 0 : ++nActive;
828 0 : }//end ipad
829 : }//ind irow
830 : }//end isec
831 0 : if (nActive>0){
832 0 : for (Int_t i=0;i<npar;++i){
833 0 : pulserQdeviations.GetMatrixArray()[i]/=nActive;
834 0 : varQMean/=nActive;
835 : }
836 0 : }
837 0 : }
838 : //_____________________________________________________________________________________
839 : void AliTPCcalibDButil::UpdatePulserOutlierMap()
840 : {
841 : /// Update the outlier map of the pulser data
842 :
843 0 : PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
844 0 : }
845 : //_____________________________________________________________________________________
846 : void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
847 : {
848 : /// Update the outlier map of the pulser reference data
849 :
850 0 : PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
851 0 : }
852 : //_____________________________________________________________________________________
853 : void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
854 : {
855 : /// Create a map that contains outliers from the Pulser calibration data.
856 : /// The outliers include masked channels, edge pads and pads with
857 : /// too large timing and charge variations.
858 : /// fNpulserOutliers is the number of outliers in the Pulser calibration data.
859 : /// those do not contain masked and edge pads
860 :
861 0 : if (!pulT||!pulQ) {
862 : //reset map
863 0 : pulOut->Multiply(0.);
864 0 : fNpulserOutliers=-1;
865 0 : return;
866 : }
867 : AliTPCCalROC *rocMasked=0x0;
868 0 : fNpulserOutliers=0;
869 :
870 : //Create Outlier Map
871 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
872 0 : AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
873 0 : AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
874 0 : AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
875 0 : if (!tmeanROC||!qmeanROC) {
876 : //reset outliers in this ROC
877 0 : outROC->Multiply(0.);
878 0 : continue;
879 : }
880 0 : if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
881 : // Double_t dummy=0;
882 : // Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
883 : // Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
884 0 : UInt_t nrows=tmeanROC->GetNrows();
885 0 : for (UInt_t irow=0;irow<nrows;++irow){
886 0 : UInt_t npads=tmeanROC->GetNPads(irow);
887 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
888 : Int_t outlier=0,masked=0;
889 0 : Float_t q=qmeanROC->GetValue(irow,ipad);
890 0 : Float_t t=tmeanROC->GetValue(irow,ipad);
891 : //masked channels are outliers
892 0 : if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
893 : //edge pads are outliers
894 0 : if (ipad==0||ipad==npads-1) masked=1;
895 : //channels with too large charge or timing deviation from the meadian are outliers
896 : // if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
897 0 : if (q<fPulQminLimit && !masked) outlier=1;
898 : //check for nan
899 0 : if ( !(q<10000000) || !(t<10000000)) outlier=1;
900 0 : outROC->SetValue(irow,ipad,outlier+masked);
901 0 : fNpulserOutliers+=outlier;
902 : }
903 : }
904 0 : }
905 0 : }
906 : //_____________________________________________________________________________________
907 : AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
908 : {
909 : /// Create pad time0 object from pulser and/or CE data, depending on the selected model
910 : /// Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
911 : /// Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
912 : /// Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
913 : ///
914 : /// In case model 2 is invoked - gy arival time gradient is also returned
915 :
916 0 : gyA=0;
917 0 : gyC=0;
918 0 : AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
919 : // decide between different models
920 0 : if (model==0||model==1){
921 0 : TVectorD vMean;
922 0 : if (model==1) ProcessPulser(vMean);
923 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
924 0 : AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
925 0 : if (!rocPulTmean) continue;
926 0 : AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
927 0 : AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
928 0 : Float_t mean=rocPulTmean->GetMean(rocOut);
929 : //treat case where a whole partition is masked
930 0 : if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
931 0 : if (model==1) {
932 0 : Int_t type=isec/18;
933 0 : mean=vMean[type];
934 0 : }
935 0 : UInt_t nrows=rocTime0->GetNrows();
936 0 : for (UInt_t irow=0;irow<nrows;++irow){
937 0 : UInt_t npads=rocTime0->GetNPads(irow);
938 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
939 0 : Float_t time=rocPulTmean->GetValue(irow,ipad);
940 : //in case of an outlier pad use the mean of the altro values.
941 : //This should be the most precise guess in that case.
942 0 : if (rocOut->GetValue(irow,ipad)) {
943 0 : time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
944 0 : if ( TMath::Abs(time)<kAlmost0 ) time=mean;
945 : }
946 0 : Float_t val=time-mean;
947 0 : rocTime0->SetValue(irow,ipad,val);
948 : }
949 : }
950 0 : }
951 0 : } else if (model==2){
952 0 : Double_t pgya,pgyc,pchi2a,pchi2c;
953 0 : AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
954 0 : fCETmean->Add(padPulser,-1.);
955 0 : TVectorD vA,vC;
956 0 : AliTPCCalPad outCE("outCE","outCE");
957 0 : Int_t nOut;
958 0 : ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
959 0 : AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
960 : // AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
961 0 : if (!padFit) { delete padPulser; return 0;}
962 0 : gyA=vA[2];
963 0 : gyC=vC[2];
964 0 : fCETmean->Add(padPulser,1.);
965 0 : padTime0->Add(fCETmean);
966 0 : padTime0->Add(padFit,-1);
967 0 : delete padPulser;
968 0 : TVectorD vFitROC;
969 0 : TMatrixD mFitROC;
970 0 : Float_t chi2;
971 0 : for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
972 0 : AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
973 0 : AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
974 0 : AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
975 0 : AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
976 0 : rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
977 0 : AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
978 0 : Float_t mean=rocPulTmean->GetMean(rocOutPul);
979 0 : if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
980 0 : UInt_t nrows=rocTime0->GetNrows();
981 0 : for (UInt_t irow=0;irow<nrows;++irow){
982 0 : UInt_t npads=rocTime0->GetNPads(irow);
983 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
984 0 : Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
985 0 : if (rocOutCE->GetValue(irow,ipad)){
986 0 : Float_t valOut=rocCEfit->GetValue(irow,ipad);
987 0 : if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
988 0 : rocTime0->SetValue(irow,ipad,valOut);
989 0 : }
990 : }
991 : }
992 0 : delete rocCEfit;
993 : }
994 0 : delete padFit;
995 0 : }
996 0 : Double_t median = padTime0->GetMedian();
997 0 : padTime0->Add(-median); // normalize to median
998 : return padTime0;
999 0 : }
1000 : //_____________________________________________________________________________________
1001 : Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *const rocOut)
1002 : {
1003 : /// GetMeanAlto information
1004 :
1005 0 : if (roc==0) return 0.;
1006 0 : const Int_t sector=roc->GetSector();
1007 0 : AliTPCROC *tpcRoc=AliTPCROC::Instance();
1008 0 : const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
1009 : Float_t mean=0;
1010 : Int_t n=0;
1011 :
1012 : //loop over a small range around the requested pad (+-10 rows/pads)
1013 0 : for (Int_t irow=row-10;irow<row+10;++irow){
1014 0 : if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
1015 0 : for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
1016 0 : if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
1017 0 : const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
1018 0 : if (altroRoc!=altroCurr) continue;
1019 0 : if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
1020 0 : Float_t val=roc->GetValue(irow,ipad);
1021 0 : mean+=val;
1022 0 : ++n;
1023 0 : }
1024 0 : }
1025 0 : if (n>0) mean/=n;
1026 : return mean;
1027 0 : }
1028 : //_____________________________________________________________________________________
1029 : void AliTPCcalibDButil::SetRefFile(const char* filename)
1030 : {
1031 : /// load cal pad objects form the reference file
1032 :
1033 0 : TDirectory *currDir=gDirectory;
1034 0 : TFile f(filename);
1035 0 : fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
1036 0 : fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
1037 : //pulser data
1038 0 : fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
1039 0 : fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
1040 0 : fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
1041 : //CE data
1042 0 : fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
1043 0 : fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
1044 0 : fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
1045 : //Altro data
1046 : // fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
1047 : // fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
1048 : // fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
1049 : // fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
1050 0 : fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
1051 0 : f.Close();
1052 0 : currDir->cd();
1053 0 : }
1054 : //_____________________________________________________________________________________
1055 : void AliTPCcalibDButil::UpdateRefDataFromOCDB()
1056 : {
1057 : /// set reference data from OCDB Reference map
1058 :
1059 0 : if (!fRefMap) {
1060 0 : AliWarning("Referenc map not set!");
1061 0 : return;
1062 : }
1063 :
1064 0 : TString cdbPath;
1065 : AliCDBEntry* entry = 0x0;
1066 : Bool_t hasAnyChanged=kFALSE;
1067 :
1068 : //pedestals
1069 0 : cdbPath="TPC/Calib/Pedestals";
1070 0 : if (HasRefChanged(cdbPath.Data())){
1071 : hasAnyChanged=kTRUE;
1072 : //delete old entries
1073 0 : if (fRefPedestals) delete fRefPedestals;
1074 0 : if (fRefPedestalMasked) delete fRefPedestalMasked;
1075 0 : fRefPedestals=fRefPedestalMasked=0x0;
1076 : //get new entries
1077 0 : entry=GetRefEntry(cdbPath.Data());
1078 0 : if (entry){
1079 0 : entry->SetOwner(kTRUE);
1080 0 : fRefPedestals=GetRefCalPad(entry);
1081 0 : delete entry;
1082 0 : fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1083 0 : }
1084 : }
1085 :
1086 : //noise
1087 0 : cdbPath="TPC/Calib/PadNoise";
1088 0 : if (HasRefChanged(cdbPath.Data())){
1089 : hasAnyChanged=kTRUE;
1090 : //delete old entry
1091 0 : if (fRefPadNoise) delete fRefPadNoise;
1092 0 : fRefPadNoise=0x0;
1093 : //get new entry
1094 0 : entry=GetRefEntry(cdbPath.Data());
1095 0 : if (entry){
1096 0 : entry->SetOwner(kTRUE);
1097 0 : fRefPadNoise=GetRefCalPad(entry);
1098 0 : delete entry;
1099 : }
1100 : }
1101 :
1102 : //pulser
1103 0 : cdbPath="TPC/Calib/Pulser";
1104 0 : if (HasRefChanged(cdbPath.Data())){
1105 : hasAnyChanged=kTRUE;
1106 : //delete old entries
1107 0 : if (fRefPulserTmean) delete fRefPulserTmean;
1108 0 : if (fRefPulserTrms) delete fRefPulserTrms;
1109 0 : if (fRefPulserQmean) delete fRefPulserQmean;
1110 0 : if (fRefPulserMasked) delete fRefPulserMasked;
1111 0 : fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1112 : //get new entries
1113 0 : entry=GetRefEntry(cdbPath.Data());
1114 0 : if (entry){
1115 0 : entry->SetOwner(kTRUE);
1116 0 : fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1117 0 : fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1118 0 : fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1119 0 : delete entry;
1120 0 : fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1121 0 : }
1122 : }
1123 :
1124 : //ce
1125 0 : cdbPath="TPC/Calib/CE";
1126 0 : if (HasRefChanged(cdbPath.Data())){
1127 : hasAnyChanged=kTRUE;
1128 : //delete old entries
1129 0 : if (fRefCETmean) delete fRefCETmean;
1130 0 : if (fRefCETrms) delete fRefCETrms;
1131 0 : if (fRefCEQmean) delete fRefCEQmean;
1132 0 : if (fRefCEMasked) delete fRefCEMasked;
1133 0 : fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1134 : //get new entries
1135 0 : entry=GetRefEntry(cdbPath.Data());
1136 0 : if (entry){
1137 0 : entry->SetOwner(kTRUE);
1138 0 : fRefCETmean=GetRefCalPad(entry,"CETmean");
1139 0 : fRefCETrms=GetRefCalPad(entry,"CETrms");
1140 0 : fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1141 0 : delete entry;
1142 0 : fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1143 0 : }
1144 : }
1145 :
1146 : //altro data
1147 0 : cdbPath="TPC/Calib/AltroConfig";
1148 0 : if (HasRefChanged(cdbPath.Data())){
1149 : hasAnyChanged=kTRUE;
1150 : //delete old entries
1151 0 : if (fRefALTROFPED) delete fRefALTROFPED;
1152 0 : if (fRefALTROZsThr) delete fRefALTROZsThr;
1153 0 : if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1154 0 : if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1155 0 : if (fRefALTROMasked) delete fRefALTROMasked;
1156 0 : fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1157 : //get new entries
1158 0 : entry=GetRefEntry(cdbPath.Data());
1159 0 : if (entry){
1160 0 : entry->SetOwner(kTRUE);
1161 0 : fRefALTROFPED=GetRefCalPad(entry,"FPED");
1162 0 : fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1163 0 : fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1164 0 : fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1165 0 : fRefALTROMasked=GetRefCalPad(entry,"Masked");
1166 0 : delete entry;
1167 : }
1168 : }
1169 :
1170 : //raw data
1171 : /*
1172 : cdbPath="TPC/Calib/Raw";
1173 : if (HasRefChanged(cdbPath.Data())){
1174 : hasAnyChanged=kTRUE;
1175 : //delete old entry
1176 : if (fRefCalibRaw) delete fRefCalibRaw;
1177 : //get new entry
1178 : entry=GetRefEntry(cdbPath.Data());
1179 : if (entry){
1180 : entry->SetOwner(kTRUE);
1181 : TObjArray *arr=(TObjArray*)entry->GetObject();
1182 : if (!arr){
1183 : AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1184 : } else {
1185 : fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1186 : }
1187 : }
1188 : }
1189 : */
1190 :
1191 : //data qa
1192 0 : cdbPath="TPC/Calib/QA";
1193 0 : if (HasRefChanged(cdbPath.Data())){
1194 : hasAnyChanged=kTRUE;
1195 : //delete old entry
1196 0 : if (fRefDataQA) delete fRefDataQA;
1197 : //get new entry
1198 0 : entry=GetRefEntry(cdbPath.Data());
1199 0 : if (entry){
1200 0 : entry->SetOwner(kTRUE);
1201 0 : fRefDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1202 0 : if (!fRefDataQA){
1203 0 : AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1204 : } else {
1205 0 : fRefDataQA=(AliTPCdataQA*)fRefDataQA->Clone();
1206 : }
1207 0 : delete entry;
1208 : }
1209 : }
1210 :
1211 :
1212 : //update current reference maps
1213 0 : if (hasAnyChanged){
1214 0 : if (fCurrentRefMap) delete fCurrentRefMap;
1215 0 : fCurrentRefMap=(TMap*)fRefMap->Clone();
1216 0 : }
1217 0 : }
1218 : //_____________________________________________________________________________________
1219 : AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1220 : {
1221 : /// TObjArray object type case
1222 : /// find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1223 :
1224 : AliTPCCalPad *pad=0x0;
1225 0 : TObjArray *arr=(TObjArray*)entry->GetObject();
1226 0 : if (!arr){
1227 0 : AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1228 0 : return pad;
1229 : }
1230 0 : pad=(AliTPCCalPad*)arr->FindObject(objName);
1231 0 : if (!pad) {
1232 0 : AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
1233 0 : return pad;
1234 : }
1235 0 : return (AliTPCCalPad*)pad->Clone();
1236 0 : }
1237 : //_____________________________________________________________________________________
1238 : AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1239 : {
1240 : /// AliTPCCalPad object type case
1241 : /// cast object to a calPad and store it in 'pad'
1242 :
1243 0 : AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1244 0 : if (!pad) {
1245 0 : AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1246 0 : return 0x0;
1247 : }
1248 0 : pad=(AliTPCCalPad*)pad->Clone();
1249 0 : return pad;
1250 0 : }
1251 : //_____________________________________________________________________________________
1252 : AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1253 : {
1254 : /// set altro masked channel map for 'cdbPath'
1255 :
1256 : AliTPCCalPad* pad=0x0;
1257 0 : const Int_t run=GetReferenceRun(cdbPath);
1258 0 : if (run<0) {
1259 0 : AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1260 0 : return pad;
1261 : }
1262 0 : AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1263 0 : if (!entry) {
1264 0 : AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1265 0 : return pad;
1266 : }
1267 0 : pad=GetRefCalPad(entry,"Masked");
1268 0 : if (pad) pad->SetNameTitle(name,name);
1269 0 : entry->SetOwner(kTRUE);
1270 0 : delete entry;
1271 0 : return pad;
1272 0 : }
1273 : //_____________________________________________________________________________________
1274 : void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1275 : /// Get Reference map
1276 :
1277 0 : if (run<0) run=fCalibDB->GetRun();
1278 0 : TString cdbPath="TPC/Calib/Ref";
1279 0 : AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1280 0 : if (!entry) {
1281 0 : AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
1282 0 : fRefMap=0;
1283 0 : return;
1284 : }
1285 0 : entry->SetOwner(kTRUE);
1286 0 : fRefMap=(TMap*)(entry->GetObject());
1287 0 : AliCDBId &id=entry->GetId();
1288 0 : fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1289 0 : }
1290 : //_____________________________________________________________________________________
1291 : Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1292 : {
1293 : /// check whether a reference cdb entry has changed
1294 :
1295 0 : if (!fCurrentRefMap) return kTRUE;
1296 0 : if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1297 0 : return kFALSE;
1298 0 : }
1299 : //_____________________________________________________________________________________
1300 : AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1301 : {
1302 : /// get the reference AliCDBEntry for 'cdbPath'
1303 :
1304 0 : const Int_t run=GetReferenceRun(cdbPath);
1305 0 : if (run<0) {
1306 0 : AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1307 0 : return 0;
1308 : }
1309 0 : AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1310 0 : if (!entry) {
1311 0 : AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1312 0 : return 0;
1313 : }
1314 0 : return entry;
1315 0 : }
1316 : //_____________________________________________________________________________________
1317 : Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
1318 : /// Get reference run number for the specified OCDB path
1319 :
1320 0 : if (!fCurrentRefMap) return -2;
1321 0 : TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1322 0 : if (!str) return -2;
1323 0 : return (Int_t)str->GetString().Atoi();
1324 0 : }
1325 : //_____________________________________________________________________________________
1326 : Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
1327 : /// Get reference run number for the specified OCDB path
1328 :
1329 0 : if (!fRefMap) return -1;
1330 0 : TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1331 0 : if (!str) return -1;
1332 0 : return (Int_t)str->GetString().Atoi();
1333 0 : }
1334 : //_____________________________________________________________________________________
1335 : AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
1336 : /// Author: marian.ivanov@cern.ch
1337 : ///
1338 : /// Create outlier map for CE study
1339 : /// Parameters:
1340 : /// Return value - outlyer map
1341 : /// noutlyersCE - number of outlyers
1342 : /// minSignal - minimal total Q signal
1343 : /// cutRMSMin - minimal width of the signal in respect to the median
1344 : /// cutRMSMax - maximal width of the signal in respect to the median
1345 : /// cutMaxDistT - maximal deviation from time median per chamber
1346 : ///
1347 : /// Outlyers criteria:
1348 : /// 0. Exclude masked pads
1349 : /// 1. Exclude first two rows in IROC and last two rows in OROC
1350 : /// 2. Exclude edge pads
1351 : /// 3. Exclude channels with too large variations
1352 : /// 4. Exclude pads with too small signal
1353 : /// 5. Exclude signal with outlyers RMS
1354 : /// 6. Exclude channels to far from the chamber median
1355 :
1356 0 : noutliersCE=0;
1357 : //create outlier map
1358 : AliTPCCalPad *out=ceOut;
1359 0 : if (!out) out= new AliTPCCalPad("outCE","outCE");
1360 : AliTPCCalROC *rocMasked=0x0;
1361 0 : if (!fCETmean) return 0;
1362 0 : if (!fCETrms) return 0;
1363 0 : if (!fCEQmean) return 0;
1364 : //
1365 : //loop over all channels
1366 : //
1367 0 : Double_t rmsMedian = fCETrms->GetMedian();
1368 0 : for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1369 0 : AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
1370 0 : if (!rocData) continue;
1371 0 : if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1372 0 : AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1373 0 : AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
1374 0 : AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1375 0 : Double_t trocMedian = rocData->GetMedian();
1376 : //
1377 0 : if (!rocData || !rocCEQ || !rocCETrms || !rocData) {
1378 0 : noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1379 0 : rocOut->Add(1.);
1380 0 : continue;
1381 : }
1382 : //
1383 : //select outliers
1384 0 : UInt_t nrows=rocData->GetNrows();
1385 0 : for (UInt_t irow=0;irow<nrows;++irow){
1386 0 : UInt_t npads=rocData->GetNPads(irow);
1387 0 : for (UInt_t ipad=0;ipad<npads;++ipad){
1388 0 : rocOut->SetValue(irow,ipad,0);
1389 0 : Float_t valTmean=rocData->GetValue(irow,ipad);
1390 0 : Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1391 0 : Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1392 : //0. exclude masked pads
1393 0 : if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1394 0 : rocOut->SetValue(irow,ipad,1);
1395 0 : continue;
1396 : }
1397 : //1. exclude first two rows in IROC and last two rows in OROC
1398 0 : if (iroc<36){
1399 0 : if (irow<2) rocOut->SetValue(irow,ipad,1);
1400 : } else {
1401 0 : if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1402 : }
1403 : //2. exclude edge pads
1404 0 : if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1405 : //exclude values that are exactly 0
1406 0 : if ( TMath::Abs(valTmean)<kAlmost0) {
1407 0 : rocOut->SetValue(irow,ipad,1);
1408 0 : ++noutliersCE;
1409 0 : }
1410 : //3. exclude channels with too large variations
1411 0 : if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1412 0 : rocOut->SetValue(irow,ipad,1);
1413 0 : ++noutliersCE;
1414 0 : }
1415 : //
1416 : //4. exclude channels with too small signal
1417 0 : if (valQmean<minSignal) {
1418 0 : rocOut->SetValue(irow,ipad,1);
1419 0 : ++noutliersCE;
1420 0 : }
1421 : //
1422 : //5. exclude channels with too small rms
1423 0 : if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1424 0 : rocOut->SetValue(irow,ipad,1);
1425 0 : ++noutliersCE;
1426 0 : }
1427 : //
1428 : //6. exclude channels to far from the chamber median
1429 0 : if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1430 0 : rocOut->SetValue(irow,ipad,1);
1431 0 : ++noutliersCE;
1432 0 : }
1433 0 : }
1434 : }
1435 0 : }
1436 : //
1437 : return out;
1438 0 : }
1439 :
1440 :
1441 : AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1442 : /// Author: marian.ivanov@cern.ch
1443 : ///
1444 : /// Create outlier map for Pulser
1445 : /// Parameters:
1446 : /// Return value - outlyer map
1447 : /// noutlyersPulser - number of outlyers
1448 : /// cutTime - absolute cut - distance to the median of chamber
1449 : /// cutnRMSQ - nsigma cut from median q distribution per chamber
1450 : /// cutnRMSrms - nsigma cut from median rms distribution
1451 : /// Outlyers criteria:
1452 : /// 0. Exclude masked pads
1453 : /// 1. Exclude time outlyers (default 3 time bins)
1454 : /// 2. Exclude q outlyers (default 5 sigma)
1455 : /// 3. Exclude rms outlyers (default 5 sigma)
1456 :
1457 0 : noutliersPulser=0;
1458 : AliTPCCalPad *out=pulserOut;
1459 0 : if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1460 : AliTPCCalROC *rocMasked=0x0;
1461 0 : if (!fPulserTmean) return 0;
1462 0 : if (!fPulserTrms) return 0;
1463 0 : if (!fPulserQmean) return 0;
1464 : //
1465 : //loop over all channels
1466 : //
1467 0 : for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1468 0 : if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1469 0 : AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1470 0 : AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1471 0 : AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1472 0 : AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1473 : //
1474 0 : Double_t rocMedianT = rocData->GetMedian();
1475 0 : Double_t rocMedianQ = rocPulserQ->GetMedian();
1476 0 : Double_t rocRMSQ = rocPulserQ->GetRMS();
1477 0 : Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1478 0 : Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1479 0 : for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1480 0 : rocOut->SetValue(ichannel,0);
1481 0 : Float_t valTmean=rocData->GetValue(ichannel);
1482 0 : Float_t valQmean=rocPulserQ->GetValue(ichannel);
1483 0 : Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1484 : Float_t valMasked =0;
1485 0 : if (rocMasked) valMasked = rocMasked->GetValue(ichannel);
1486 : Int_t isOut=0;
1487 0 : if (valMasked>0.5) isOut=1;
1488 0 : if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1489 0 : if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1490 0 : if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1491 0 : rocOut->SetValue(ichannel,isOut);
1492 0 : if (isOut) noutliersPulser++;
1493 : }
1494 : }
1495 0 : return out;
1496 0 : }
1497 :
1498 :
1499 : AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1500 : /// Author : Marian Ivanov
1501 : /// Create pad time0 correction map using information from the CE and from pulser
1502 : ///
1503 : /// Return PadTime0 to be used for time0 relative alignment
1504 : /// if dump file specified intermediat results are dumped to the fiel and can be visualized
1505 : /// using $ALICE_ROOT/TPC/script/gui application
1506 : ///
1507 : /// fitResultsA - fitParameters A side
1508 : /// fitResultsC - fitParameters C side
1509 : /// chi2A - chi2/ndf for A side (assuming error 1 time bin)
1510 : /// chi2C - chi2/ndf for C side (assuming error 1 time bin)
1511 : ///
1512 : /// Algorithm:
1513 : /// 1. Find outlier map for CE
1514 : /// 2. Find outlier map for Pulser
1515 : /// 3. Replace outlier by median at given sector (median without outliers)
1516 : /// 4. Substract from the CE data pulser
1517 : /// 5. Fit the CE with formula
1518 : /// 5.1) (IROC-OROC) offset
1519 : /// 5.2) gx
1520 : /// 5.3) gy
1521 : /// 5.4) (lx-xmid)
1522 : /// 5.5) (IROC-OROC)*(lx-xmid)
1523 : /// 5.6) (ly/lx)^2
1524 : /// 6. Substract gy fit dependence from the CE data
1525 : /// 7. Add pulser back to CE data
1526 : /// 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1527 : /// 9. return CE data
1528 : ///
1529 : /// Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1530 : /// Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1531 :
1532 : // fit formula
1533 : const char *formulaIn="(-1.+2.*(sector<36))*0.5++gx++gy++(lx-134.)++(-1.+2.*(sector<36))*0.5*(lx-134)++((ly/lx)^2/(0.1763)^2)";
1534 : // output for fit formula
1535 : const char *formulaAll="1++(-1.+2.*(sector<36))*0.5++gx++gy++(lx-134.)++(-1.+2.*(sector<36))*0.5*(lx-134)++((ly/lx)^2/(0.1763)^2)";
1536 : // gy part of formula
1537 : const char *formulaOut="0++0*(-1.+2.*(sector<36))*0.5++0*gx++gy++0*(lx-134.)++0*(-1.+2.*(sector<36))*0.5*(lx-134)++0*((ly/lx)^2/(0.1763)^2)";
1538 : //
1539 : //
1540 0 : if (!fCETmean) return 0;
1541 0 : Double_t pgya,pgyc,pchi2a,pchi2c;
1542 0 : AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1543 0 : AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1544 :
1545 0 : AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1546 0 : AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1547 0 : AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1548 0 : AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1549 0 : padPulser->SetName("padPulser");
1550 0 : padPulserOut->SetName("padPulserOut");
1551 0 : padCE->SetName("padCE");
1552 0 : padCEIn->SetName("padCEIn");
1553 0 : padCEOut->SetName("padCEOut");
1554 0 : padOut->SetName("padOut");
1555 :
1556 : //
1557 : // make combined outlyers map
1558 : // and replace outlyers in maps with median for chamber
1559 : //
1560 0 : for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1561 0 : AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1562 0 : AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1563 0 : AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1564 0 : AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1565 0 : AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1566 0 : Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1567 0 : Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1568 0 : for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1569 0 : if (rocPulserOut->GetValue(ichannel)>0) {
1570 0 : rocPulser->SetValue(ichannel,pulserMedian);
1571 0 : rocOut->SetValue(ichannel,1);
1572 0 : }
1573 0 : if (rocCEOut->GetValue(ichannel)>0) {
1574 0 : rocCE->SetValue(ichannel,ceMedian);
1575 0 : rocOut->SetValue(ichannel,1);
1576 0 : }
1577 : }
1578 : }
1579 : //
1580 : // remove pulser time 0
1581 : //
1582 0 : padCE->Add(padPulser,-1);
1583 : //
1584 : // Make fits
1585 : //
1586 0 : TMatrixD dummy;
1587 0 : Float_t chi2Af,chi2Cf;
1588 0 : padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1589 0 : chi2A=chi2Af;
1590 0 : chi2C=chi2Cf;
1591 : //
1592 0 : AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1593 0 : padCEFitGY->SetName("padCEFitGy");
1594 : //
1595 0 : AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1596 0 : padCEFit->SetName("padCEFit");
1597 : //
1598 0 : AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1599 0 : padCEDiff->SetName("padCEDiff");
1600 0 : padCEDiff->Add(padCEFit,-1.);
1601 : //
1602 : //
1603 0 : padCE->Add(padCEFitGY,-1.);
1604 :
1605 0 : padCE->Add(padPulser,1.);
1606 0 : Double_t padmedian = padCE->GetMedian();
1607 0 : padCE->Add(-padmedian); // normalize to median
1608 : //
1609 : // Replace outliers by fit value - median of diff per given chamber -GY fit
1610 : //
1611 0 : for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1612 0 : AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1613 0 : AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1614 0 : AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1615 0 : AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1616 0 : AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1617 : //
1618 0 : Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1619 0 : for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1620 0 : if (rocOut->GetValue(ichannel)==0) continue;
1621 0 : Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1622 0 : rocCE->SetValue(ichannel,value);
1623 0 : }
1624 : }
1625 : //
1626 : //
1627 0 : if (dumpfile){
1628 : //dump to the file - result can be visualized
1629 0 : AliTPCPreprocessorOnline preprocesor;
1630 0 : preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1631 0 : preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1632 0 : preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1633 0 : preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1634 : //
1635 0 : preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1636 0 : preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1637 : //
1638 0 : preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1639 0 : preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1640 0 : preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1641 0 : preprocesor.DumpToFile(dumpfile);
1642 0 : }
1643 0 : delete padPulser;
1644 0 : delete padPulserOut;
1645 0 : delete padCEIn;
1646 0 : delete padCEOut;
1647 0 : delete padOut;
1648 0 : delete padCEDiff;
1649 0 : delete padCEFitGY;
1650 : return padCE;
1651 0 : }
1652 :
1653 :
1654 :
1655 :
1656 :
1657 : Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1658 : /// find the closest point to xref in x direction
1659 : /// return dx and value
1660 :
1661 0 : dx = 0;
1662 0 : y = 0;
1663 :
1664 0 : if(!graph) return 0;
1665 0 : if(graph->GetN() < 1) return 0;
1666 :
1667 : Int_t index=0;
1668 0 : index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1669 0 : if (index<0) index=0;
1670 0 : if(graph->GetN()==1) {
1671 0 : dx = xref-graph->GetX()[index];
1672 0 : }
1673 : else {
1674 0 : if (index>=graph->GetN()-1) index=graph->GetN()-2;
1675 0 : if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1676 0 : dx = xref-graph->GetX()[index];
1677 : }
1678 0 : y = graph->GetY()[index];
1679 : return index;
1680 0 : }
1681 :
1682 : Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1683 : /// Get the correction of the trigger offset
1684 : /// combining information from the laser track calibration
1685 : /// and from cosmic calibration
1686 : ///
1687 : /// run - run number
1688 : /// timeStamp - tim stamp in seconds
1689 : /// deltaT - integration period to calculate offset
1690 : /// deltaTLaser -max validity of laser data
1691 : /// valType - 0 - median, 1- mean
1692 : ///
1693 : /// Integration vaues are just recomendation - if not possible to get points
1694 : /// automatically increase the validity by factor 2
1695 : /// (recursive algorithm until one month of data taking)
1696 :
1697 : const Float_t kLaserCut=0.0005;
1698 : const Int_t kMaxPeriod=3600*24*30*12; // one year max
1699 : const Int_t kMinPoints=20;
1700 : //
1701 0 : TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1702 0 : if (!array) {
1703 0 : AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1704 0 : }
1705 0 : array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1706 0 : if (!array) return 0;
1707 : //
1708 : TGraphErrors *laserA[3]={0,0,0};
1709 : TGraphErrors *laserC[3]={0,0,0};
1710 : TGraphErrors *cosmicAll=0;
1711 0 : laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1712 0 : laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1713 0 : cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1714 : //
1715 : //
1716 0 : if (!cosmicAll) return 0;
1717 0 : Int_t nmeasC=cosmicAll->GetN();
1718 0 : Float_t *tdelta = new Float_t[nmeasC];
1719 : Int_t nused=0;
1720 0 : for (Int_t i=0;i<nmeasC;i++){
1721 0 : if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1722 0 : Float_t ccosmic=cosmicAll->GetY()[i];
1723 0 : Double_t yA=0,yC=0,dA=0,dC=0;
1724 0 : if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1725 0 : if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1726 : //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1727 : //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1728 : //
1729 0 : if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1730 : Float_t claser=0;
1731 0 : if (TMath::Abs(yA-yC)<kLaserCut) {
1732 0 : claser=(yA-yC)*0.5;
1733 0 : }else{
1734 0 : if (i%2==0) claser=yA;
1735 0 : if (i%2==1) claser=yC;
1736 : }
1737 0 : tdelta[nused]=ccosmic-claser;
1738 0 : nused++;
1739 0 : }
1740 0 : if (nused<kMinPoints &&deltaT<kMaxPeriod) {
1741 0 : delete [] tdelta;
1742 0 : return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1743 : }
1744 0 : if (nused<kMinPoints) {
1745 0 : delete [] tdelta;
1746 : //AliWarning("AliFatal: No time offset calibration available\n");
1747 0 : return 0;
1748 : }
1749 0 : Double_t median = TMath::Median(nused,tdelta);
1750 0 : Double_t mean = TMath::Mean(nused,tdelta);
1751 0 : delete [] tdelta;
1752 0 : return (valType==0) ? median:mean;
1753 0 : }
1754 :
1755 : Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1756 : /// Get the correction of the drift velocity
1757 : /// combining information from the laser track calibration
1758 : /// and from cosmic calibration
1759 : ///
1760 : /// dist - return value - distance to closest point in graph
1761 : /// run - run number
1762 : /// timeStamp - tim stamp in seconds
1763 : /// deltaT - integration period to calculate time0 offset
1764 : /// deltaTLaser -max validity of laser data
1765 : /// valType - 0 - median, 1- mean
1766 : ///
1767 : /// Integration vaues are just recomendation - if not possible to get points
1768 : /// automatically increase the validity by factor 2
1769 : /// (recursive algorithm until one month of data taking)
1770 :
1771 4 : TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1772 2 : if (!array) {
1773 2 : AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1774 2 : }
1775 2 : array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1776 4 : if (!array) return 0;
1777 : TGraphErrors *cosmicAll=0;
1778 0 : cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1779 0 : if (!cosmicAll) return 0;
1780 0 : Double_t grY=0;
1781 0 : AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1782 :
1783 0 : Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1784 0 : Double_t vcosmic = AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1785 0 : if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1786 0 : if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1787 0 : return vcosmic-t0;
1788 :
1789 : /*
1790 : Example usage:
1791 :
1792 : Int_t run=89000
1793 : TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1794 : cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1795 : laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1796 : //
1797 : Double_t *yvd= new Double_t[cosmicAll->GetN()];
1798 : Double_t *yt0= new Double_t[cosmicAll->GetN()];
1799 : for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1800 : for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1801 :
1802 : TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1803 : TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1804 :
1805 : */
1806 :
1807 2 : }
1808 :
1809 : const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1810 : {
1811 : /// Create a default name for the gui file
1812 :
1813 0 : return Form("guiRefTreeRun%s.root",GetRefValidity());
1814 : }
1815 :
1816 : Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1817 : {
1818 : /// Create a gui reference tree
1819 : /// if dirname and filename are empty default values will be used
1820 : /// this is the recommended way of using this function
1821 : /// it allows to check whether a file with the given run validity alredy exists
1822 :
1823 0 : if (!AliCDBManager::Instance()->GetDefaultStorage()){
1824 0 : AliError("Default Storage not set. Cannot create reference calibration Tree!");
1825 0 : return kFALSE;
1826 : }
1827 :
1828 0 : TString file=filename;
1829 0 : if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1830 :
1831 0 : AliTPCPreprocessorOnline prep;
1832 : //noise and pedestals
1833 0 : if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1834 0 : if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1835 0 : if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1836 : //pulser data
1837 0 : if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1838 0 : if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1839 0 : if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1840 0 : if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1841 : //CE data
1842 0 : if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1843 0 : if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1844 0 : if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1845 0 : if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1846 : //Altro data
1847 0 : if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1848 0 : if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1849 0 : if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1850 0 : if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1851 0 : if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1852 : //QA
1853 0 : AliTPCdataQA *dataQA=fRefDataQA;
1854 0 : if (dataQA) {
1855 0 : if (dataQA->GetNLocalMaxima())
1856 0 : prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1857 0 : if (dataQA->GetMaxCharge())
1858 0 : prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1859 0 : if (dataQA->GetMeanCharge())
1860 0 : prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1861 0 : if (dataQA->GetNoThreshold())
1862 0 : prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1863 0 : if (dataQA->GetNTimeBins())
1864 0 : prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1865 0 : if (dataQA->GetNPads())
1866 0 : prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1867 0 : if (dataQA->GetTimePosition())
1868 0 : prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1869 : }
1870 0 : prep.DumpToFile(file.Data());
1871 : return kTRUE;
1872 0 : }
1873 :
1874 : Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1875 : /// Get the correction of the drift velocity using the offline laser tracks calbration
1876 : ///
1877 : /// run - run number
1878 : /// timeStamp - tim stamp in seconds
1879 : /// deltaT - integration period to calculate time0 offset
1880 : /// side - 0 - A side, 1 - C side, 2 - mean from both sides
1881 : /// Note in case no data form both A and C side - the value from active side used
1882 :
1883 4 : TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1884 :
1885 2 : return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1886 : }
1887 :
1888 : Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1889 : /// Get the correction of the drift velocity using the online laser tracks calbration
1890 : ///
1891 : /// run - run number
1892 : /// timeStamp - tim stamp in seconds
1893 : /// deltaT - integration period to calculate time0 offset
1894 : /// side - 0 - A side, 1 - C side, 2 - mean from both sides
1895 : /// Note in case no data form both A and C side - the value from active side used
1896 :
1897 0 : TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1898 :
1899 0 : Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1900 0 : AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1901 0 : if (!param) return 0;
1902 :
1903 : //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1904 0 : dv*=param->GetDriftV()/2.61301900000000000e+06;
1905 0 : if (dv>1e-20) dv=1/dv-1;
1906 0 : else return 0;
1907 : // T/P correction
1908 0 : TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1909 :
1910 0 : AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
1911 0 : AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
1912 :
1913 : Double_t corrPTA=0;
1914 : Double_t corrPTC=0;
1915 :
1916 0 : if (temp&&press) {
1917 0 : AliTPCCalibVdrift corr(temp,press,0);
1918 0 : corrPTA=corr.GetPTRelative(timeStamp,0);
1919 0 : corrPTC=corr.GetPTRelative(timeStamp,1);
1920 0 : }
1921 :
1922 0 : if (side==0) dv -= corrPTA;
1923 0 : if (side==1) dv -= corrPTC;
1924 0 : if (side==2) dv -= (corrPTA+corrPTC)/2;
1925 :
1926 : return dv;
1927 0 : }
1928 :
1929 : Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1930 : Int_t side, TObjArray * const array){
1931 : /// common drift velocity retrieval for online and offline method
1932 :
1933 : TGraphErrors *grlaserA=0;
1934 : TGraphErrors *grlaserC=0;
1935 : Double_t vlaserA=0, vlaserC=0;
1936 6 : if (!array) return 0;
1937 0 : grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1938 0 : grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1939 0 : Double_t deltaY;
1940 0 : if (grlaserA && grlaserA->GetN()>0) {
1941 0 : AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1942 0 : if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1943 0 : else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1944 : }
1945 0 : if (grlaserC && grlaserC->GetN()>0) {
1946 0 : AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1947 0 : if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1948 0 : else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1949 : }
1950 0 : if (side==0) return vlaserA;
1951 0 : if (side==1) return vlaserC;
1952 0 : Double_t mdrift=(vlaserA+vlaserC)*0.5;
1953 0 : if (!grlaserA) return vlaserC;
1954 0 : if (!grlaserC) return vlaserA;
1955 0 : return mdrift;
1956 2 : }
1957 :
1958 :
1959 : Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1960 : /// Get the correction of the drift velocity using the CE laser data
1961 : /// combining information from the CE, laser track calibration
1962 : /// and P/T calibration
1963 : ///
1964 : /// run - run number
1965 : /// timeStamp - tim stamp in seconds
1966 : /// deltaT - integration period to calculate time0 offset
1967 : /// side - 0 - A side, 1 - C side, 2 - mean from both sides
1968 :
1969 4 : TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
1970 2 : if (!arrT) return 0;
1971 2 : AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1972 2 : TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1973 2 : AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1974 : //
1975 : //
1976 : Double_t corrPTA = 0, corrPTC=0;
1977 : Double_t ltime0A = 0, ltime0C=0;
1978 2 : Double_t gry=0;
1979 : Double_t corrA=0, corrC=0;
1980 : Double_t timeA=0, timeC=0;
1981 : const Double_t kEpsilon = 0.00001;
1982 2 : TGraph *graphA = (TGraph*)arrT->At(72);
1983 2 : TGraph *graphC = (TGraph*)arrT->At(73);
1984 2 : if (!graphA && !graphC) return 0.;
1985 4 : if (graphA &&graphA->GetN()>0) {
1986 0 : AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1987 0 : timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1988 0 : Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1989 0 : ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1990 0 : if(ltime0A < kEpsilon) return 0;
1991 0 : if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1992 0 : corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1993 0 : corrA-=corrPTA;
1994 0 : }
1995 4 : if (graphC&&graphC->GetN()>0){
1996 0 : AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1997 0 : timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1998 0 : Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1999 0 : ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
2000 0 : if(ltime0C < kEpsilon) return 0;
2001 0 : if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
2002 0 : corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
2003 0 : corrC-=corrPTC;
2004 0 : }
2005 :
2006 2 : if (side ==0 ) return corrA;
2007 2 : if (side ==1 ) return corrC;
2008 2 : Double_t corrM= (corrA+corrC)*0.5;
2009 2 : if (!graphA) corrM=corrC;
2010 2 : if (!graphC) corrM=corrA;
2011 : return corrM;
2012 4 : }
2013 :
2014 : Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2015 : /// return drift velocity using the TPC-ITS matchin method
2016 : /// return also distance to the closest point
2017 :
2018 4 : TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2019 : TGraphErrors *graph=0;
2020 2 : dist=0;
2021 4 : if (!array) return 0;
2022 : //array->ls();
2023 0 : graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
2024 0 : if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_TOFB_TPC_DRIFTVD");
2025 0 : if (!graph) return 0;
2026 0 : Double_t deltaY;
2027 0 : AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2028 0 : Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2029 : return value;
2030 2 : }
2031 :
2032 : Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2033 : /// Get time dependent time 0 (trigger delay in cm) correction
2034 : /// Arguments:
2035 : /// timestamp - timestamp
2036 : /// run - run number
2037 : ///
2038 : /// Notice - Extrapolation outside of calibration range - using constant function
2039 :
2040 4 : TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2041 : TGraphErrors *graph=0;
2042 2 : dist=0;
2043 4 : if (!array) return 0;
2044 0 : graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
2045 0 : if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
2046 0 : if (!graph) return 0;
2047 0 : Double_t deltaY;
2048 0 : AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2049 0 : Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2050 : return value;
2051 2 : }
2052 :
2053 :
2054 :
2055 :
2056 :
2057 : Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2058 : /// VERY obscure method - we need something in framework
2059 : /// Find the TPC runs with temperature OCDB entry
2060 : /// cache the start and end of the run
2061 :
2062 0 : AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2063 0 : if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2064 0 : if (!storage) return 0;
2065 0 : TString path=storage->GetURI();
2066 0 : TString runsT;
2067 : {
2068 0 : TString command;
2069 0 : if (path.Contains("local")){ // find the list if local system
2070 0 : path.ReplaceAll("local://","");
2071 0 : path+="TPC/Calib/Temperature";
2072 0 : command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2073 : }
2074 0 : runsT=gSystem->GetFromPipe(command);
2075 0 : }
2076 0 : TObjArray *arr= runsT.Tokenize("r");
2077 0 : if (!arr) return 0;
2078 : //
2079 0 : TArrayI indexes(arr->GetEntries());
2080 0 : TArrayI runs(arr->GetEntries());
2081 : Int_t naccept=0;
2082 0 : {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2083 0 : Int_t irunN = atoi(arr->At(irun)->GetName());
2084 0 : if (irunN<startRun) continue;
2085 0 : if (irunN>stopRun) continue;
2086 0 : runs[naccept]=irunN;
2087 0 : naccept++;
2088 0 : }}
2089 0 : delete arr;
2090 0 : fRuns.Set(naccept);
2091 0 : fRunsStart.Set(fRuns.fN);
2092 0 : fRunsStop.Set(fRuns.fN);
2093 0 : TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2094 0 : for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2095 :
2096 : //
2097 : AliCDBEntry * entry = 0;
2098 0 : {for (Int_t irun=0;irun<fRuns.fN; irun++){
2099 0 : entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2100 0 : if (!entry) continue;
2101 0 : AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2102 0 : if (!tmpRun) continue;
2103 0 : fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2104 0 : fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2105 : //AliInfo(Form("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec()));
2106 0 : }}
2107 : return fRuns.fN;
2108 0 : }
2109 :
2110 :
2111 : Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2112 : /// binary search - find the run for given time stamp
2113 :
2114 0 : Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2115 0 : Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2116 : Int_t cindex = -1;
2117 0 : for (Int_t index=index0; index<=index1; index++){
2118 0 : if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2119 0 : if (debug) {
2120 0 : AliInfo(Form("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime));
2121 0 : }
2122 : }
2123 0 : if (cindex<0) cindex =(index0+index1)/2;
2124 0 : if (cindex<0) {
2125 0 : return 0;
2126 : }
2127 0 : return fRuns[cindex];
2128 0 : }
2129 :
2130 :
2131 :
2132 :
2133 :
2134 : TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2135 : /// filter outlyer measurement
2136 : /// Only points around median +- sigmaCut filtered
2137 :
2138 0 : if (!graph) return 0;
2139 : Int_t kMinPoints=2;
2140 0 : Int_t npoints0 = graph->GetN();
2141 : Int_t npoints=0;
2142 : Float_t rmsY=0;
2143 : //
2144 : //
2145 0 : if (npoints0<kMinPoints) return 0;
2146 :
2147 0 : Double_t *outx=new Double_t[npoints0];
2148 0 : Double_t *outy=new Double_t[npoints0];
2149 0 : for (Int_t iter=0; iter<3; iter++){
2150 : npoints=0;
2151 0 : for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2152 0 : if (graph->GetY()[ipoint]==0) continue;
2153 0 : if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2154 0 : outx[npoints] = graph->GetX()[ipoint];
2155 0 : outy[npoints] = graph->GetY()[ipoint];
2156 0 : npoints++;
2157 0 : }
2158 0 : if (npoints<=1) break;
2159 0 : medianY =TMath::Median(npoints,outy);
2160 0 : rmsY =TMath::RMS(npoints,outy);
2161 : }
2162 : TGraph *graphOut=0;
2163 0 : if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2164 0 : delete [] outx;
2165 0 : delete [] outy;
2166 : return graphOut;
2167 0 : }
2168 :
2169 :
2170 : TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2171 : /// filter outlyer measurement
2172 : /// Only points around median +- cut filtered
2173 :
2174 0 : if (!graph) return 0;
2175 : Int_t kMinPoints=2;
2176 0 : Int_t npoints0 = graph->GetN();
2177 : Int_t npoints=0;
2178 : // Float_t rmsY=0;
2179 : //
2180 : //
2181 0 : if (npoints0<kMinPoints) return 0;
2182 :
2183 0 : Double_t *outx=new Double_t[npoints0];
2184 0 : Double_t *outy=new Double_t[npoints0];
2185 0 : for (Int_t iter=0; iter<3; iter++){
2186 : npoints=0;
2187 0 : for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2188 0 : if (graph->GetY()[ipoint]==0) continue;
2189 0 : if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2190 0 : outx[npoints] = graph->GetX()[ipoint];
2191 0 : outy[npoints] = graph->GetY()[ipoint];
2192 0 : npoints++;
2193 0 : }
2194 0 : if (npoints<=1) break;
2195 0 : medianY =TMath::Median(npoints,outy);
2196 : // rmsY =TMath::RMS(npoints,outy);
2197 : }
2198 : TGraph *graphOut=0;
2199 0 : if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2200 0 : delete [] outx;
2201 0 : delete [] outy;
2202 : return graphOut;
2203 0 : }
2204 :
2205 :
2206 :
2207 : TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
2208 : /// filter outlyer measurement
2209 : /// Only points with normalized errors median +- sigmaCut filtered
2210 :
2211 : Int_t kMinPoints=10;
2212 0 : Int_t npoints0 = graph->GetN();
2213 : Int_t npoints=0;
2214 : Float_t medianErr=0, rmsErr=0;
2215 : //
2216 : //
2217 0 : if (npoints0<kMinPoints) return 0;
2218 :
2219 0 : Double_t *outx=new Double_t[npoints0];
2220 0 : Double_t *outy=new Double_t[npoints0];
2221 0 : Double_t *erry=new Double_t[npoints0];
2222 0 : Double_t *nerry=new Double_t[npoints0];
2223 0 : Double_t *errx=new Double_t[npoints0];
2224 :
2225 0 : for (Int_t iter=0; iter<3; iter++){
2226 : npoints=0;
2227 0 : for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2228 0 : nerry[npoints] = graph->GetErrorY(ipoint);
2229 0 : if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2230 0 : erry[npoints] = graph->GetErrorY(ipoint);
2231 0 : outx[npoints] = graph->GetX()[ipoint];
2232 0 : outy[npoints] = graph->GetY()[ipoint];
2233 0 : errx[npoints] = graph->GetErrorY(ipoint);
2234 0 : npoints++;
2235 0 : }
2236 0 : if (npoints==0) break;
2237 0 : medianErr=TMath::Median(npoints,erry);
2238 0 : medianY =TMath::Median(npoints,outy);
2239 0 : rmsErr =TMath::RMS(npoints,erry);
2240 : }
2241 : TGraphErrors *graphOut=0;
2242 0 : if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2243 0 : delete []outx;
2244 0 : delete []outy;
2245 0 : delete []erry;
2246 0 : delete []nerry;
2247 0 : delete []errx;
2248 : return graphOut;
2249 0 : }
2250 :
2251 :
2252 : void AliTPCcalibDButil::Sort(TGraph *graph){
2253 : /// sort array - neccessay for approx
2254 :
2255 0 : Int_t npoints = graph->GetN();
2256 0 : Int_t *indexes=new Int_t[npoints];
2257 0 : Double_t *outx=new Double_t[npoints];
2258 0 : Double_t *outy=new Double_t[npoints];
2259 0 : TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2260 0 : for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2261 0 : for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2262 0 : for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2263 0 : for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2264 :
2265 0 : delete [] indexes;
2266 0 : delete [] outx;
2267 0 : delete [] outy;
2268 0 : }
2269 : void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2270 : /// smmoth graph - mean on the interval
2271 :
2272 0 : Sort(graph);
2273 0 : Int_t npoints = graph->GetN();
2274 0 : Double_t *outy=new Double_t[npoints];
2275 :
2276 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2277 0 : Double_t lx=graph->GetX()[ipoint];
2278 0 : Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2279 0 : Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2280 0 : if (index0<0) index0=0;
2281 0 : if (index1>=npoints-1) index1=npoints-1;
2282 0 : if ((index1-index0)>1){
2283 0 : outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2284 0 : }else{
2285 0 : outy[ipoint]=graph->GetY()[ipoint];
2286 : }
2287 : }
2288 : // TLinearFitter fitter(3,"pol2");
2289 : // for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2290 : // Double_t lx=graph->GetX()[ipoint];
2291 : // Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2292 : // Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2293 : // if (index0<0) index0=0;
2294 : // if (index1>=npoints-1) index1=npoints-1;
2295 : // fitter.ClearPoints();
2296 : // for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2297 : // if ((index1-index0)>1){
2298 : // outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2299 : // }else{
2300 : // outy[ipoint]=graph->GetY()[ipoint];
2301 : // }
2302 : // }
2303 :
2304 :
2305 :
2306 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2307 0 : graph->GetY()[ipoint] = outy[ipoint];
2308 : }
2309 0 : delete[] outy;
2310 0 : }
2311 :
2312 : Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
2313 : /// Use constant interpolation outside of range
2314 :
2315 0 : if (!graph) {
2316 0 : AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2317 0 : return 0;
2318 : }
2319 :
2320 0 : if (graph->GetN()<1){
2321 0 : AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
2322 0 : return 0;
2323 : }
2324 :
2325 :
2326 0 : if (xref<graph->GetX()[0]) return graph->GetY()[0];
2327 0 : if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2328 :
2329 : // AliInfo(Form("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref)));
2330 :
2331 0 : if(graph->GetN()==1)
2332 0 : return graph->Eval(graph->GetX()[0]);
2333 :
2334 :
2335 0 : return graph->Eval(xref);
2336 0 : }
2337 :
2338 : Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2339 : /// Use constant interpolation outside of range also for spline fits
2340 :
2341 0 : if (!graph) {
2342 0 : AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2343 0 : return 0;
2344 : }
2345 0 : if (graph->GetKnots()<1){
2346 0 : AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph");
2347 0 : return 0;
2348 : }
2349 0 : if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2350 0 : if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2351 0 : return graph->Eval( xref);
2352 0 : }
2353 :
2354 : Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2355 : /// Filter DCS sensor information
2356 : /// ymin - minimal value
2357 : /// ymax - max value
2358 : /// maxdy - maximal deirivative
2359 : /// sigmaCut - cut on values and derivative in terms of RMS distribution
2360 : /// Return value - accepted fraction
2361 : ///
2362 : /// Algorithm:
2363 : ///
2364 : /// 0. Calculate median and rms of values in specified range
2365 : /// 1. Filter out outliers - median+-sigmaCut*rms
2366 : /// values replaced by median
2367 :
2368 0 : AliSplineFit * fit = sensor->GetFit();
2369 0 : if (!fit) return 0.;
2370 0 : Int_t nknots = fit->GetKnots();
2371 0 : if (nknots==0) {
2372 0 : delete fit;
2373 0 : sensor->SetFit(0);
2374 0 : return 0;
2375 : }
2376 : //
2377 0 : Double_t *yin0 = new Double_t[nknots];
2378 0 : Double_t *yin1 = new Double_t[nknots];
2379 : Int_t naccept=0;
2380 :
2381 0 : for (Int_t iknot=0; iknot< nknots; iknot++){
2382 0 : if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2383 0 : yin0[naccept] = fit->GetY0()[iknot];
2384 0 : yin1[naccept] = fit->GetY1()[iknot];
2385 0 : if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2386 0 : naccept++;
2387 0 : }
2388 : }
2389 0 : if (naccept<1) {
2390 0 : delete fit;
2391 0 : sensor->SetFit(0);
2392 0 : delete [] yin0;
2393 0 : delete [] yin1;
2394 0 : return 0.;
2395 : }
2396 :
2397 : Double_t medianY0=0, medianY1=0;
2398 : Double_t rmsY0 =0, rmsY1=0;
2399 0 : medianY0 = TMath::Median(naccept, yin0);
2400 0 : medianY1 = TMath::Median(naccept, yin1);
2401 0 : rmsY0 = TMath::RMS(naccept, yin0);
2402 0 : rmsY1 = TMath::RMS(naccept, yin1);
2403 : naccept=0;
2404 : //
2405 : // 1. Filter out outliers - median+-sigmaCut*rms
2406 : // values replaced by median
2407 : // if replaced the derivative set to 0
2408 : //
2409 0 : for (Int_t iknot=0; iknot< nknots; iknot++){
2410 : Bool_t isOK=kTRUE;
2411 0 : if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2412 0 : if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2413 0 : if (nknots<2) fit->GetY1()[iknot]=0;
2414 0 : if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2415 0 : if (!isOK){
2416 0 : fit->GetY0()[iknot]=medianY0;
2417 0 : fit->GetY1()[iknot]=0;
2418 0 : }else{
2419 0 : naccept++;
2420 : }
2421 : }
2422 0 : delete [] yin0;
2423 0 : delete [] yin1;
2424 0 : return Float_t(naccept)/Float_t(nknots);
2425 0 : }
2426 :
2427 : Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2428 : /// Filter temperature array
2429 : /// tempArray - array of temperatures -
2430 : /// ymin - minimal accepted temperature - default 15
2431 : /// ymax - maximal accepted temperature - default 22
2432 : /// sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2433 : /// return value - fraction of filtered sensors
2434 :
2435 : const Double_t kMaxDy=0.1;
2436 0 : Int_t nsensors=tempArray->NumSensors();
2437 0 : if (nsensors==0) return 0.;
2438 : Int_t naccept=0;
2439 0 : for (Int_t isensor=0; isensor<nsensors; isensor++){
2440 0 : AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2441 0 : if (!sensor) continue;
2442 0 : FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2443 0 : if (sensor->GetFit()==0){
2444 : //delete sensor;
2445 0 : tempArray->RemoveSensorNum(isensor);
2446 0 : }else{
2447 0 : naccept++;
2448 : }
2449 0 : }
2450 0 : return Float_t(naccept)/Float_t(nsensors);
2451 0 : }
2452 :
2453 :
2454 : void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
2455 : /// Filter CE data
2456 : /// Input parameters:
2457 : /// deltaT - smoothing window (in seconds)
2458 : /// cutAbs - max distance of the time info to the median (in time bins)
2459 : /// cutSigma - max distance (in the RMS)
2460 : /// pcstream - optional debug streamer to store original and filtered info
2461 : /// Hardwired parameters:
2462 : /// kMinPoints =10; // minimal number of points to define the CE
2463 : /// kMinSectors=12; // minimal number of sectors to define sideCE
2464 : /// Algorithm:
2465 : /// 0. Filter almost emty graphs (kMinPoints=10)
2466 : /// 1. calculate median and RMS per side
2467 : /// 2. Filter graphs - in respect with side medians
2468 : /// - cutAbs and cutDelta used
2469 : /// 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2470 : /// 4. Calculate mean for A side and C side
2471 :
2472 : const Int_t kMinPoints =10; // minimal number of points to define the CE
2473 : const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2474 : const Int_t kMinTime =400; // minimal arrival time of CE
2475 0 : TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2476 0 : Double_t medianY=0;
2477 0 : TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2478 0 : if (!cearray) return;
2479 : Double_t tmin=-1;
2480 : Double_t tmax=-1;
2481 : //
2482 : //
2483 0 : AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2484 0 : AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
2485 0 : if ( tempMapCE && cavernPressureCE){
2486 : //
2487 : // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2488 : // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2489 : // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2490 : Bool_t isOK=kTRUE;
2491 0 : if (isOK) {
2492 : // recalculate P/T correction map for time of the CE
2493 0 : AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2494 0 : driftCalib->SetName("driftPTCE");
2495 0 : driftCalib->SetTitle("driftPTCE");
2496 0 : cearray->AddLast(driftCalib);
2497 0 : }
2498 0 : }
2499 : //
2500 : // 0. Filter almost emty graphs
2501 : //
2502 :
2503 0 : for (Int_t i=0; i<72;i++){
2504 0 : TGraph *graph= (TGraph*)arrT->At(i);
2505 0 : if (!graph) continue;
2506 0 : graph->Sort();
2507 0 : if (graph->GetN()<kMinPoints){
2508 0 : arrT->AddAt(0,i);
2509 0 : delete graph; // delete empty graph
2510 0 : continue;
2511 : }
2512 0 : if (tmin<0) tmin = graph->GetX()[0];
2513 0 : if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2514 : //
2515 0 : if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2516 0 : if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2517 0 : }
2518 : //
2519 : // 1. calculate median and RMS per side
2520 : //
2521 0 : TArrayF arrA(100000), arrC(100000);
2522 0 : Int_t nA=0, nC=0;
2523 : Double_t medianA=0, medianC=0;
2524 : Double_t rmsA=0, rmsC=0;
2525 0 : for (Int_t isec=0; isec<72;isec++){
2526 0 : TGraph *graph= (TGraph*)arrT->At(isec);
2527 0 : if (!graph) continue;
2528 0 : for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2529 0 : if (graph->GetY()[ipoint]<kMinTime) continue;
2530 0 : if (nA>=arrA.fN) arrA.Set(nA*2);
2531 0 : if (nC>=arrC.fN) arrC.Set(nC*2);
2532 0 : if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2533 0 : if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2534 : }
2535 0 : }
2536 0 : if (nA>0){
2537 0 : medianA=TMath::Median(nA,arrA.fArray);
2538 0 : rmsA =TMath::RMS(nA,arrA.fArray);
2539 0 : }
2540 0 : if (nC>0){
2541 0 : medianC=TMath::Median(nC,arrC.fArray);
2542 0 : rmsC =TMath::RMS(nC,arrC.fArray);
2543 0 : }
2544 : //
2545 : // 2. Filter graphs - in respect with side medians
2546 : //
2547 0 : TArrayD vecX(100000), vecY(100000);
2548 0 : for (Int_t isec=0; isec<72;isec++){
2549 0 : TGraph *graph= (TGraph*)arrT->At(isec);
2550 0 : if (!graph) continue;
2551 0 : Double_t median = (isec%36<18) ? medianA: medianC;
2552 0 : Double_t rms = (isec%36<18) ? rmsA: rmsC;
2553 : Int_t naccept=0;
2554 : // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2555 0 : for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
2556 0 : if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2557 0 : if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2558 0 : vecX[naccept]= graph->GetX()[ipoint];
2559 0 : vecY[naccept]= graph->GetY()[ipoint];
2560 0 : naccept++;
2561 0 : }
2562 0 : if (naccept<kMinPoints){
2563 0 : arrT->AddAt(0,isec);
2564 0 : delete graph; // delete empty graph
2565 0 : continue;
2566 : }
2567 0 : TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2568 0 : delete graph;
2569 0 : arrT->AddAt(graph2,isec);
2570 0 : }
2571 : //
2572 : // 3. Cut in respect wit the graph median
2573 : //
2574 0 : for (Int_t i=0; i<72;i++){
2575 0 : TGraph *graph= (TGraph*)arrT->At(i);
2576 0 : if (!graph) continue;
2577 : //
2578 : // filter in range
2579 : //
2580 0 : TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2581 0 : if (!graphTS0) continue;
2582 0 : if (graphTS0->GetN()<kMinPoints) {
2583 0 : delete graphTS0;
2584 0 : delete graph;
2585 0 : arrT->AddAt(0,i);
2586 0 : continue;
2587 : }
2588 0 : TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2589 0 : if (!graphTS) continue;
2590 0 : graphTS->Sort();
2591 0 : AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2592 0 : if (pcstream){
2593 0 : Int_t run = AliTPCcalibDB::Instance()->GetRun();
2594 0 : (*pcstream)<<"filterCE"<<
2595 0 : "run="<<run<<
2596 0 : "isec="<<i<<
2597 0 : "mY="<<medianY<<
2598 0 : "graph.="<<graph<<
2599 0 : "graphTS0.="<<graphTS0<<
2600 0 : "graphTS.="<<graphTS<<
2601 : "\n";
2602 0 : }
2603 0 : delete graphTS0;
2604 0 : arrT->AddAt(graphTS,i);
2605 0 : delete graph;
2606 0 : }
2607 : //
2608 : // Recalculate the mean time A side C side
2609 : //
2610 0 : TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2611 0 : Int_t meanPoints=(nA+nC)/72; // mean number of points
2612 0 : for (Int_t itime=0; itime<200; itime++){
2613 0 : nA=0, nC=0;
2614 0 : Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2615 0 : for (Int_t i=0; i<72;i++){
2616 0 : TGraph *graph= (TGraph*)arrT->At(i);
2617 0 : if (!graph) continue;
2618 0 : if (graph->GetN()<(meanPoints/4)) continue;
2619 0 : if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2620 0 : if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2621 0 : }
2622 0 : xA[itime]=time;
2623 0 : xC[itime]=time;
2624 0 : yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2625 0 : yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2626 0 : eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2627 0 : eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2628 : }
2629 : //
2630 0 : Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2631 0 : Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2632 0 : if (pcstream){
2633 0 : Int_t run = AliTPCcalibDB::Instance()->GetRun();
2634 0 : (*pcstream)<<"filterAC"<<
2635 0 : "run="<<run<<
2636 0 : "nA="<<nA<<
2637 0 : "nC="<<nC<<
2638 0 : "rmsTA="<<rmsTA<<
2639 0 : "rmsTC="<<rmsTC<<
2640 : "\n";
2641 0 : }
2642 : //
2643 0 : TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2644 0 : TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2645 0 : TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2646 0 : if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2647 0 : TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2648 0 : if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2649 0 : delete grA;
2650 0 : delete grC;
2651 0 : if (nA<kMinSectors) arrT->AddAt(0,72);
2652 0 : else arrT->AddAt(graphTSA,72);
2653 0 : if (nC<kMinSectors) arrT->AddAt(0,73);
2654 0 : else arrT->AddAt(graphTSC,73);
2655 0 : }
2656 :
2657 :
2658 : void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
2659 : /// Filter Drift velocity measurement using the tracks
2660 : /// 0. remove outlyers - error based
2661 : /// cutSigma
2662 :
2663 : const Int_t kMinPoints=1; // minimal number of points to define value
2664 0 : TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2665 0 : Double_t medianY=0;
2666 0 : if (!arrT) return;
2667 0 : for (Int_t i=0; i<arrT->GetEntries();i++){
2668 0 : TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
2669 0 : if (!graph) continue;
2670 0 : if (graph->GetN()<kMinPoints){
2671 0 : delete graph;
2672 0 : arrT->AddAt(0,i);
2673 0 : continue;
2674 : }
2675 : TGraphErrors *graph2 = NULL;
2676 0 : if(graph->GetN()<10) {
2677 0 : graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2678 0 : if (!graph2) {
2679 0 : delete graph; arrT->AddAt(0,i); continue;
2680 : }
2681 : }
2682 : else {
2683 0 : graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2684 0 : if (!graph2) {
2685 0 : delete graph; arrT->AddAt(0,i); continue;
2686 : }
2687 : }
2688 0 : if (graph2->GetN()<1) {
2689 0 : delete graph; arrT->AddAt(0,i); continue;
2690 : }
2691 0 : graph2->SetName(graph->GetName());
2692 0 : graph2->SetTitle(graph->GetTitle());
2693 0 : arrT->AddAt(graph2,i);
2694 0 : if (pcstream){
2695 0 : (*pcstream)<<"filterTracks"<<
2696 0 : "run="<<run<<
2697 0 : "isec="<<i<<
2698 0 : "mY="<<medianY<<
2699 0 : "graph.="<<graph<<
2700 0 : "graph2.="<<graph2<<
2701 : "\n";
2702 0 : }
2703 0 : delete graph;
2704 0 : }
2705 0 : }
2706 :
2707 :
2708 :
2709 :
2710 :
2711 : Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2712 : /// get laser time offset
2713 : /// median around timeStamp+-deltaT
2714 : /// QA - chi2 needed for later usage - to be added
2715 : /// - currently cut on error
2716 :
2717 : Int_t kMinPoints=1;
2718 : Double_t kMinDelay=0.01;
2719 : Double_t kMinDelayErr=0.0001;
2720 :
2721 0 : TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2722 0 : if (!array) return 0;
2723 : TGraphErrors *tlaser=0;
2724 0 : if (array){
2725 0 : if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2726 0 : if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2727 : }
2728 0 : if (!tlaser) return 0;
2729 0 : Int_t npoints0= tlaser->GetN();
2730 0 : if (npoints0==0) return 0;
2731 0 : Double_t *xlaser = new Double_t[npoints0];
2732 0 : Double_t *ylaser = new Double_t[npoints0];
2733 : Int_t npoints=0;
2734 0 : for (Int_t i=0;i<npoints0;i++){
2735 : //printf("%d\n",i);
2736 0 : if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2737 0 : if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2738 0 : xlaser[npoints]=tlaser->GetX()[npoints];
2739 0 : ylaser[npoints]=tlaser->GetY()[npoints];
2740 0 : npoints++;
2741 0 : }
2742 : //
2743 : //
2744 0 : Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2745 0 : Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2746 : //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2747 0 : if (index0<0) index0=0;
2748 0 : if (index1>=npoints-1) index1=npoints-1;
2749 0 : if (index1-index0<kMinPoints) {
2750 0 : delete [] ylaser;
2751 0 : delete [] xlaser;
2752 0 : return 0;
2753 : }
2754 : //
2755 : //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2756 0 : Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2757 0 : delete [] ylaser;
2758 0 : delete [] xlaser;
2759 : return mean;
2760 0 : }
2761 :
2762 :
2763 :
2764 :
2765 : void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
2766 : /// Filter Goofie data
2767 : /// goofieArray - points will be filtered
2768 : /// deltaT - smmothing time window
2769 : /// cutSigma - outler sigma cut in rms
2770 : /// minVn, maxVd- range absolute cut for variable vd/pt
2771 : /// - to be tuned
2772 : ///
2773 : /// Ignore goofie if not enough points
2774 :
2775 : const Int_t kMinPoints = 3;
2776 : //
2777 :
2778 0 : TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2779 0 : TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2780 0 : TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2781 0 : TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2782 0 : if (!graphvd) return;
2783 0 : if (graphvd->GetN()<kMinPoints){
2784 0 : delete graphvd;
2785 0 : goofieArray->GetSensorNum(2)->SetGraph(0);
2786 0 : return;
2787 : }
2788 : //
2789 : // 1. Caluclate medians of critical variables
2790 : // drift velcocity
2791 : // P/T
2792 : // area near peak
2793 : // area far peak
2794 : //
2795 0 : Double_t medianpt=0;
2796 0 : Double_t medianvd=0, sigmavd=0;
2797 0 : Double_t medianan=0;
2798 0 : Double_t medianaf=0;
2799 0 : Int_t entries=graphvd->GetN();
2800 0 : Double_t yvdn[10000];
2801 : Int_t nvd=0;
2802 : //
2803 0 : for (Int_t ipoint=0; ipoint<entries; ipoint++){
2804 0 : if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2805 0 : if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2806 0 : if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2807 0 : yvdn[nvd++]=graphvd->GetY()[ipoint];
2808 0 : }
2809 0 : if (nvd<kMinPoints){
2810 0 : delete graphvd;
2811 0 : goofieArray->GetSensorNum(2)->SetGraph(0);
2812 0 : return;
2813 : }
2814 : //
2815 0 : Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2816 0 : if (nuni>=kMinPoints){
2817 0 : AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2818 0 : }else{
2819 0 : medianvd = TMath::Median(nvd, yvdn);
2820 : }
2821 :
2822 0 : TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2823 0 : TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2824 0 : TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2825 0 : TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2826 0 : TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2827 0 : TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2828 0 : delete graphpt0;
2829 0 : delete graphpt1;
2830 0 : delete graphan0;
2831 0 : delete graphan1;
2832 0 : delete graphaf0;
2833 0 : delete graphaf1;
2834 : //
2835 : // 2. Make outlyer graph
2836 : //
2837 : Int_t nOK=0;
2838 0 : TGraph graphOut(*graphvd);
2839 0 : for (Int_t i=0; i<entries;i++){
2840 : //
2841 : Bool_t isOut=kFALSE;
2842 0 : if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2843 0 : if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2844 :
2845 0 : if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2846 0 : isOut|=kTRUE;
2847 0 : if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2848 0 : if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2849 0 : if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2850 0 : graphOut.GetY()[i]= (isOut)?1:0;
2851 0 : if (!isOut) nOK++;
2852 0 : }
2853 0 : if (nOK<kMinPoints) {
2854 0 : delete graphvd;
2855 0 : goofieArray->GetSensorNum(2)->SetGraph(0);
2856 0 : return;
2857 : }
2858 : //
2859 : // 3. Filter out outlyers - and smooth
2860 : //
2861 0 : TVectorF vmedianArray(goofieArray->NumSensors());
2862 0 : TVectorF vrmsArray(goofieArray->NumSensors());
2863 0 : Double_t xnew[10000];
2864 0 : Double_t ynew[10000];
2865 0 : TObjArray junk;
2866 0 : junk.SetOwner(kTRUE);
2867 0 : Bool_t isOK=kTRUE;
2868 : //
2869 : //
2870 0 : for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2871 0 : isOK=kTRUE;
2872 0 : AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2873 : TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2874 : //
2875 0 : if (!sensor) continue;
2876 0 : graphOld = sensor->GetGraph();
2877 0 : if (graphOld) {
2878 0 : sensor->SetGraph(0);
2879 : Int_t nused=0;
2880 0 : for (Int_t i=0;i<entries;i++){
2881 0 : if (graphOut.GetY()[i]>0.5) continue;
2882 0 : xnew[nused]=graphOld->GetX()[i];
2883 0 : ynew[nused]=graphOld->GetY()[i];
2884 0 : nused++;
2885 0 : }
2886 0 : graphNew = new TGraph(nused,xnew,ynew);
2887 0 : junk.AddLast(graphNew);
2888 0 : junk.AddLast(graphOld);
2889 0 : Double_t median=0;
2890 0 : graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2891 0 : if (graphNew0!=0){
2892 0 : junk.AddLast(graphNew0);
2893 0 : graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2894 0 : if (graphNew1!=0){
2895 0 : junk.AddLast(graphNew1);
2896 0 : graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2897 0 : if (graphNew2!=0) {
2898 0 : vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2899 0 : AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2900 0 : AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2901 0 : AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2902 : // AliInfo(Form("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]));
2903 0 : vmedianArray[isensor]=median;
2904 : //
2905 0 : }
2906 : }
2907 : }
2908 0 : }
2909 0 : if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2910 0 : if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2911 0 : if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2912 0 : if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2913 0 : (*pcstream)<<"goofieA"<<
2914 0 : Form("isOK_%d.=",isensor)<<isOK<<
2915 0 : Form("s_%d.=",isensor)<<sensor<<
2916 0 : Form("gr_%d.=",isensor)<<graphOld<<
2917 0 : Form("gr0_%d.=",isensor)<<graphNew0<<
2918 0 : Form("gr1_%d.=",isensor)<<graphNew1<<
2919 0 : Form("gr2_%d.=",isensor)<<graphNew2;
2920 0 : if (isOK) sensor->SetGraph(graphNew2);
2921 0 : }
2922 0 : (*pcstream)<<"goofieA"<<
2923 0 : "vmed.="<<&vmedianArray<<
2924 0 : "vrms.="<<&vrmsArray<<
2925 : "\n";
2926 0 : junk.Delete(); // delete temoprary graphs
2927 :
2928 0 : }
2929 :
2930 :
2931 :
2932 :
2933 :
2934 : TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
2935 : /// Make a statistic matrix
2936 : /// Input parameters:
2937 : /// array - TObjArray of AliRelKalmanAlign
2938 : /// minFraction - minimal ration of accepted tracks
2939 : /// minStat - minimal statistic (number of accepted tracks)
2940 : /// maxvd - maximal deviation for the 1
2941 : /// Output matrix:
2942 : /// columns - Mean, Median, RMS
2943 : /// row - parameter type (rotation[3], translation[3], drift[3])
2944 :
2945 0 : if (!array) return 0;
2946 0 : if (array->GetEntries()<=0) return 0;
2947 : // Int_t entries = array->GetEntries();
2948 0 : Int_t entriesFast = array->GetEntriesFast();
2949 0 : TVectorD state(9);
2950 0 : TVectorD *valArray[9];
2951 0 : for (Int_t i=0; i<9; i++){
2952 0 : valArray[i] = new TVectorD(entriesFast);
2953 : }
2954 : Int_t naccept=0;
2955 0 : for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
2956 0 : AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
2957 0 : if (!kalman) continue;
2958 0 : if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
2959 0 : if (kalman->GetNUpdates()<minStat) continue;
2960 0 : if (Float_t(kalman->GetNUpdates())/Float_t(kalman->GetNTracks())<minFraction) continue;
2961 0 : kalman->GetState(state);
2962 0 : for (Int_t ipar=0; ipar<9; ipar++)
2963 0 : (*valArray[ipar])[naccept]=state[ipar];
2964 0 : naccept++;
2965 0 : }
2966 : //if (naccept<2) return 0;
2967 0 : if (naccept<1) return 0;
2968 0 : TMatrixD *pstat=new TMatrixD(9,3);
2969 : TMatrixD &stat=*pstat;
2970 0 : for (Int_t ipar=0; ipar<9; ipar++){
2971 0 : stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
2972 0 : stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
2973 0 : stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
2974 : }
2975 0 : for (Int_t i=0; i<9; i++){
2976 0 : delete valArray[i];
2977 : }
2978 : return pstat;
2979 0 : }
2980 :
2981 :
2982 : TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD* statP, Bool_t direction, Float_t sigmaCut){
2983 : /// Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
2984 : /// Input:
2985 : /// array - input array
2986 : /// stat - mean parameters statistic
2987 : /// direction -
2988 : /// sigmaCut - maximal allowed deviation from mean in terms of RMS
2989 :
2990 0 : if (!array) return 0;
2991 0 : if (array->GetEntries()<=0) return 0;
2992 0 : if (!statP) return 0;
2993 : const TMatrixD& stat = *statP;
2994 : // error increase in 1 hour
2995 : const Double_t kerrsTime[9]={
2996 : 0.00001, 0.00001, 0.00001,
2997 : 0.001, 0.001, 0.001,
2998 : 0.002, 0.01, 0.001};
2999 : //
3000 : //
3001 0 : Int_t entries = array->GetEntriesFast();
3002 0 : TObjArray *sArray= new TObjArray(entries);
3003 : AliRelAlignerKalman * sKalman =0;
3004 0 : TVectorD state(9);
3005 0 : for (Int_t i=0; i<entries; i++){
3006 0 : Int_t index=(direction)? entries-i-1:i;
3007 0 : AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3008 0 : if (!kalman) continue;
3009 : Bool_t isOK=kTRUE;
3010 0 : kalman->GetState(state);
3011 0 : for (Int_t ipar=0; ipar<9; ipar++){
3012 0 : if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3013 : }
3014 0 : if (!sKalman &&isOK) {
3015 0 : sKalman=new AliRelAlignerKalman(*kalman);
3016 0 : sKalman->SetRejectOutliers(kFALSE);
3017 0 : sKalman->SetRunNumber(kalman->GetRunNumber());
3018 0 : sKalman->SetTimeStamp(kalman->GetTimeStamp());
3019 0 : }
3020 0 : if (!sKalman) continue;
3021 0 : Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
3022 0 : for (Int_t ipar=0; ipar<9; ipar++){
3023 : // (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3024 : // (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3025 : // (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3026 0 : (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3027 : }
3028 0 : sKalman->SetRunNumber(kalman->GetRunNumber());
3029 0 : if (!isOK) sKalman->SetRunNumber(0);
3030 0 : sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3031 0 : if (!isOK) continue;
3032 0 : sKalman->SetRejectOutliers(kFALSE);
3033 0 : sKalman->SetRunNumber(kalman->GetRunNumber());
3034 0 : sKalman->SetTimeStamp(kalman->GetTimeStamp());
3035 0 : sKalman->Merge(kalman);
3036 0 : sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3037 : //sKalman->Print();
3038 0 : }
3039 : return sArray;
3040 0 : }
3041 :
3042 : TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
3043 : /// Merge 2 RelKalman arrays
3044 : /// Input:
3045 : /// arrayP - rel kalman in direction plus
3046 : /// arrayM - rel kalman in direction minus
3047 :
3048 0 : if (!arrayP) return 0;
3049 0 : if (arrayP->GetEntries()<=0) return 0;
3050 0 : if (!arrayM) return 0;
3051 0 : if (arrayM->GetEntries()<=0) return 0;
3052 :
3053 0 : Int_t entries = arrayP->GetEntriesFast();
3054 0 : TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3055 :
3056 0 : for (Int_t i=0; i<entries; i++){
3057 0 : AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3058 0 : AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3059 0 : if (!kalmanP) continue;
3060 0 : if (!kalmanM) continue;
3061 :
3062 : AliRelAlignerKalman *kalman = NULL;
3063 0 : if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3064 0 : kalman = new AliRelAlignerKalman(*kalmanP);
3065 0 : kalman->Merge(kalmanM);
3066 0 : }
3067 0 : else if (kalmanP->GetRunNumber() == 0) {
3068 0 : kalman = new AliRelAlignerKalman(*kalmanM);
3069 0 : }
3070 0 : else if (kalmanM->GetRunNumber() == 0) {
3071 0 : kalman = new AliRelAlignerKalman(*kalmanP);
3072 : }
3073 : else
3074 0 : continue;
3075 :
3076 0 : array->AddAt(kalman,i);
3077 0 : }
3078 :
3079 : return array;
3080 0 : }
3081 :
3082 : //_____________________________________________________________________________________
3083 : TTree* AliTPCcalibDButil::ConnectGainTrees(TString baseDir)
3084 : {
3085 : /// baseDir: Base directory with the raw Kr calibration trees
3086 : /// and the trees from the calibQA
3087 : /// it assumes to following structure below:
3088 : /// KryptonCalib/<year>/calibKr/calibKr.<year>.<id>.root
3089 : /// calibQAdEdx/<year>/calibQA.<year>.<perid>.tree.root
3090 : /// map/treeMapping.root
3091 :
3092 :
3093 : // === add main tree, which will be a mapping file ================
3094 0 : TFile *fin = TFile::Open(Form("%s/map/treeMapping.root",baseDir.Data()));
3095 0 : gROOT->cd();
3096 0 : TTree *tMain = (TTree*)fin->Get("calPads");
3097 0 : tMain->SetAlias("Pad0","MapPadLength.fElements==7.5"); // pad types
3098 0 : tMain->SetAlias("Pad1","MapPadLength.fElements==10.0");
3099 0 : tMain->SetAlias("Pad2","MapPadLength.fElements==15.0");
3100 : //
3101 0 : tMain->SetAlias("dRowNorm0","(row.fElements/32-1)"); // normalized distance to the center of the chamber in lx (-1,1)
3102 0 : tMain->SetAlias("dRowNorm1","(row.fElements/32-1)"); //
3103 0 : tMain->SetAlias("dRowNorm2","((row.fElements-63)/16-1)"); //
3104 0 : tMain->SetAlias("dRowNorm","(row.fElements<63)*(row.fElements/32-1)+(row.fElements>63)*((row.fElements-63)/16-1)");
3105 0 : tMain->SetAlias("dPhiNorm","(ly.fElements/lx.fElements)/(pi/18)"); // normalized distance to the center of the chamber in phi (-1,1)
3106 : //
3107 0 : tMain->SetAlias("fsector","(sector%36)+(0.5+9*atan2(ly.fElements,lx.fElements)/pi)"); // float sector number
3108 0 : tMain->SetAlias("posEdge","((pi/18.)-abs(atan(ly.fElements/lx.fElements)))*lx.fElements"); //distance to the edge
3109 :
3110 : // === add the krypton calibration trees ==========================
3111 0 : TString inputTreesKrCalib = gSystem->GetFromPipe(Form("ls %s/KryptonCalib/20*/calibKr/*.tree.root",baseDir.Data()));
3112 0 : TObjArray *arrInputTreesKrCalib = inputTreesKrCalib.Tokenize("\n");
3113 :
3114 0 : for (Int_t itree=0; itree<arrInputTreesKrCalib->GetEntriesFast(); ++itree) {
3115 0 : TFile *fin2 = TFile::Open(arrInputTreesKrCalib->At(itree)->GetName());
3116 0 : TTree *tin = (TTree*)fin2->Get("calPads");
3117 0 : gROOT->cd();
3118 0 : TString friendName=gSystem->BaseName(arrInputTreesKrCalib->At(itree)->GetName());
3119 0 : friendName.ReplaceAll("calibKr.","");
3120 0 : friendName.ReplaceAll(".tree.root","");
3121 0 : friendName="Kr."+friendName;
3122 0 : tMain->AddFriend(tin,friendName.Data());
3123 :
3124 : // set aliases
3125 :
3126 : // TODO: finish implementation of alias via lists
3127 : // const Int_t nbranchAlias = 2;
3128 : // const char* branchNames[nbranchAlias]={"spectrMean","fitMean"};
3129 : // const Int_t nbranchStat = 2;
3130 : // const char* statNames[nbranchStat] = {"Median","LTM"};
3131 : //
3132 : // for (Int_t iname=0; iname<nbranchAlias; ++iname) {
3133 : // TString branchName = TString::Format("%s.%s", friendName.Data(), bNames[i]);
3134 : //
3135 : // for (Int_t istat=0; istat<nbranchStat; ++istat) {
3136 : //
3137 : // }
3138 : // }
3139 :
3140 0 : tMain->SetAlias((friendName+".spectrMean_LTMRatio").Data(),
3141 0 : TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_LTM+0))",
3142 0 : friendName.Data(),friendName.Data()).Data());
3143 :
3144 0 : tMain->SetAlias((friendName+".spectrMean_MedianRatio").Data(),
3145 0 : TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_Median+0))",
3146 0 : friendName.Data(),friendName.Data()).Data());
3147 0 : tMain->SetAlias((friendName+".spectrMean_MeanRatio").Data(),
3148 0 : TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_Mean+0))",
3149 0 : friendName.Data(),friendName.Data()).Data());
3150 0 : tMain->SetAlias((friendName+".spectrMean_MeanRatio").Data(),
3151 0 : TString::Format("(%s.spectrMean.fElements/%s.spectrMean_Mean)",
3152 0 : friendName.Data(),friendName.Data()).Data());
3153 :
3154 0 : tMain->SetAlias((friendName+".fitMean_LTMRatio").Data(),
3155 0 : TString::Format("(%s.fitMean.fElements/(%s.fitMean_LTM+0))",
3156 0 : friendName.Data(),friendName.Data()).Data());
3157 0 : tMain->SetAlias((friendName+".fitRMS_LTMRatio").Data(),
3158 0 : TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_LTM+0))",
3159 0 : friendName.Data(),friendName.Data()).Data());
3160 :
3161 0 : tMain->SetAlias((friendName+".fitMean_MedianRatio").Data(),
3162 0 : TString::Format("(%s.fitMean.fElements/(%s.fitMean_Median+0))",
3163 0 : friendName.Data(),friendName.Data()).Data());
3164 0 : tMain->SetAlias((friendName+".fitRMS_MedianRatio").Data(),
3165 0 : TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_Median+0))",
3166 0 : friendName.Data(),friendName.Data()).Data());
3167 0 : tMain->SetAlias((friendName+".fitMean_MeanRatio").Data(),
3168 0 : TString::Format("(%s.fitMean.fElements/(%s.fitMean_Mean+0))",
3169 0 : friendName.Data(),friendName.Data()).Data());
3170 0 : tMain->SetAlias((friendName+".fitRMS_MeanRatio").Data(),
3171 0 : TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_Mean+0))",
3172 0 : friendName.Data(),friendName.Data()).Data());
3173 0 : tMain->SetAlias((friendName+".fitMean_MeanRatio").Data(),
3174 0 : TString::Format("(%s.fitMean.fElements/%s.fitMean_Mean)",
3175 0 : friendName.Data(),friendName.Data()).Data());
3176 :
3177 0 : }
3178 :
3179 : // === add the calibQA trees ======================================
3180 0 : TString inputTreesQACalib = gSystem->GetFromPipe(Form("ls %s/calibQAdEdx/20*/*.tree.root",baseDir.Data()));
3181 0 : TObjArray *arrInputTreesQACalib = inputTreesQACalib.Tokenize("\n");
3182 :
3183 0 : for (Int_t itree=0; itree<arrInputTreesQACalib->GetEntriesFast(); ++itree) {
3184 0 : TFile *fin2 = TFile::Open(arrInputTreesQACalib->At(itree)->GetName());
3185 0 : TTree *tin = (TTree*)fin2->Get("calPads");
3186 0 : gROOT->cd();
3187 0 : TString friendName=gSystem->BaseName(arrInputTreesQACalib->At(itree)->GetName());
3188 0 : friendName.ReplaceAll("calibQA.","");
3189 0 : friendName.ReplaceAll(".tree.root","");
3190 0 : friendName="QA."+friendName;
3191 0 : tMain->AddFriend(tin,friendName.Data());
3192 :
3193 : // set aliases
3194 0 : tMain->SetAlias((friendName+".MaxCharge_LTMRatio").Data(),
3195 0 : TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_LTM)",
3196 0 : friendName.Data(),friendName.Data()).Data());
3197 :
3198 0 : tMain->SetAlias((friendName+".MaxCharge_MedianRatio").Data(),
3199 0 : TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_Median)",
3200 0 : friendName.Data(),friendName.Data()).Data());
3201 0 : tMain->SetAlias((friendName+".MaxCharge_MeanRatio").Data(),
3202 0 : TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_Mean)",
3203 0 : friendName.Data(),friendName.Data()).Data());
3204 :
3205 0 : tMain->SetAlias((friendName+".MeanCharge_LTMRatio").Data(),
3206 0 : TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_LTM)",
3207 0 : friendName.Data(),friendName.Data()).Data());
3208 :
3209 0 : tMain->SetAlias((friendName+".MeanCharge_MedianRatio").Data(),
3210 0 : TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_Median)",
3211 0 : friendName.Data(),friendName.Data()).Data());
3212 0 : tMain->SetAlias((friendName+".MeanCharge_MeanRatio").Data(),
3213 0 : TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_Mean)",
3214 0 : friendName.Data(),friendName.Data()).Data());
3215 :
3216 0 : }
3217 :
3218 : return tMain;
3219 0 : }
3220 :
3221 :
3222 : //_____________________________________________________________________________________
3223 : TTree* AliTPCcalibDButil::ConnectPulserTrees(TString baseDir, TTree *tMain)
3224 : {
3225 : /// baseDir: Base directory with Pulser information
3226 : /// TTrees are added to the base tree as a friend tree
3227 : ///
3228 : /// === add the calibPulser trees ======================================
3229 :
3230 0 : TString inputTreesPulserCalib = gSystem->GetFromPipe(Form("ls %s/calibPulser/20*/*.tree.root",baseDir.Data()));
3231 0 : TObjArray *arrInputTreesPulserCalib = inputTreesPulserCalib.Tokenize("\n");
3232 0 : for (Int_t itree=0; itree<arrInputTreesPulserCalib->GetEntriesFast(); ++itree) {
3233 0 : TFile *fin2 = TFile::Open(arrInputTreesPulserCalib->At(itree)->GetName());
3234 0 : TTree *tin = (TTree*)fin2->Get("calPads");
3235 0 : gROOT->cd();
3236 0 : TString friendName=gSystem->BaseName(arrInputTreesPulserCalib->At(itree)->GetName());
3237 0 : friendName.ReplaceAll("calibPulser.","");
3238 0 : friendName.ReplaceAll(".tree.root","");
3239 0 : friendName="Pulser."+friendName;
3240 0 : tMain->AddFriend(tin,friendName.Data());
3241 : // set aliases
3242 :
3243 0 : tMain->SetAlias((friendName+".CEQmean_LTMRatio").Data(),
3244 0 : TString::Format("(%s.CEQmean.fElements/%s.CEQmean_LTM)",
3245 0 : friendName.Data(),friendName.Data()).Data());
3246 0 : tMain->SetAlias((friendName+".CEQmean_MedianRatio").Data(),
3247 0 : TString::Format("(%s.CEQmean.fElements/%s.CEQmean_Median)",
3248 0 : friendName.Data(),friendName.Data()).Data());
3249 0 : tMain->SetAlias((friendName+".CEQmean_MeanRatio").Data(),
3250 0 : TString::Format("(%s.CEQmean.fElements/%s.CEQmean_Mean)",
3251 0 : friendName.Data(),friendName.Data()).Data());
3252 : //
3253 0 : tMain->SetAlias((friendName+".CETmean_LTMDelta").Data(),
3254 0 : TString::Format("(%s.CETmean.fElements-%s.CETmean_LTM)",
3255 0 : friendName.Data(),friendName.Data()).Data());
3256 0 : tMain->SetAlias((friendName+".CETmean_MedianDelta").Data(),
3257 0 : TString::Format("(%s.CETmean.fElements-%s.CETmean_Median)",
3258 0 : friendName.Data(),friendName.Data()).Data());
3259 0 : tMain->SetAlias((friendName+".CETmean_MeanDelta").Data(),
3260 0 : TString::Format("(%s.CETmean.fElements-%s.CETmean_Mean)",
3261 0 : friendName.Data(),friendName.Data()).Data());
3262 0 : }
3263 : return tMain;
3264 0 : }
3265 :
3266 :
3267 : TTree* AliTPCcalibDButil::ConnectDistortionTrees(TString baseDir, TString selection, TTree *tMain){
3268 : /// baseDir: Base directory with Distortion information
3269 : /// TTrees are added to the base tree as a friend tree
3270 : /// If base tree not provide - first tree from list is used as base
3271 : ///
3272 : /// === add the calibDistortion trees ======================================
3273 : /// TString inputTreesDistortionCalib = gSystem->GetFromPipe(Form("ls %s/calibDistortion/20*/*.tree.root",baseDir.Data()));
3274 : /// TString baseDir="$NOTES/reconstruction/distortionFit/"; TTree *tMain=0;
3275 : /// AliTPCcalibDButil::ConnectDistortionTrees("$NOTES/reconstruction/distortionFit/", "calibTimeResHisto.root", 0);
3276 :
3277 0 : TString inputTreesDistortionCalib = "";
3278 0 : if (selection.Contains(".list")){
3279 0 : inputTreesDistortionCalib=gSystem->GetFromPipe(Form("cat %s",selection.Data()));
3280 0 : }else{
3281 0 : inputTreesDistortionCalib=gSystem->GetFromPipe(Form("find %s -iname \"%s\"",baseDir.Data(),selection.Data()));
3282 : }
3283 0 : TObjArray *arrInputTreesDistortionCalib = inputTreesDistortionCalib.Tokenize("\n");
3284 : //
3285 0 : TString treePrefix="";
3286 0 : for (Int_t itree=0; itree<arrInputTreesDistortionCalib->GetEntriesFast(); ++itree) {
3287 0 : TString fstring = arrInputTreesDistortionCalib->At(itree)->GetName();
3288 0 : if (fstring.Index("Prefix:")==0){ //get tree descriptor
3289 0 : treePrefix=TString(&(fstring[7]));
3290 0 : continue;
3291 : }
3292 0 : if (fstring.Index("#")==0){ // ignore comment field
3293 0 : continue;
3294 : }
3295 0 : TFile *finput= TFile::Open(arrInputTreesDistortionCalib->At(itree)->GetName());
3296 0 : TString strFile=arrInputTreesDistortionCalib->At(itree)->GetName();
3297 0 : TObjArray *path=strFile.Tokenize("/");
3298 0 : Int_t plength=path->GetEntries();
3299 0 : if (!finput) continue;
3300 0 : TList* list = finput->GetListOfKeys();
3301 0 : Int_t nkeys=list->GetEntries();
3302 0 : for (Int_t ikey=0; ikey<nkeys; ikey++){
3303 0 : TKey * key = (TKey*)list->At(ikey);
3304 0 : if (strstr(key->GetClassName(),"TTree")==0) continue;
3305 0 : TTree * tree = dynamic_cast<TTree*>(finput->Get(list->At(ikey)->GetName()));
3306 0 : if (!tree) continue;
3307 0 : TString friendName="";
3308 0 : if (treePrefix.Length()==0) {
3309 0 : friendName=TString::Format("%s.%s.%s",path->At(plength-3)->GetName(),path->At(plength-2)->GetName(), tree->GetName());
3310 0 : }else{
3311 0 : friendName=TString::Format("%s.%s",treePrefix.Data(), tree->GetName());
3312 : }
3313 0 : ::Info("AliTPCcalibDButil::ConnectDistortionTrees","%s",friendName.Data());
3314 0 : if (tMain==0) {
3315 : tMain=tree;
3316 0 : tree = dynamic_cast<TTree*>(finput->Get(list->At(ikey)->GetName()));
3317 0 : }
3318 0 : tMain->AddFriend(tree,friendName.Data());
3319 0 : }
3320 0 : }
3321 : // tMain->SetAlias("");
3322 : return tMain;
3323 0 : }
3324 :
3325 :
3326 : //_____________________________________________________________________________________
3327 : TTree* AliTPCcalibDButil::ConnectCalPadTrees(TString baseDir, TString pattern, TTree *tMain, Bool_t checkAliases)
3328 : {
3329 : /// baseDir: Base directory with per Pad information
3330 : /// TTrees are added to the base tree as a friend tree
3331 : /// Example usage
3332 : /// TString baseDir="/hera/alice/fsozzi/summarymaps/calib2/"; // prefix directory with calibration with slash at the end
3333 : /// TString pattern="QA/*/*root";
3334 : /// TTree * tree = AliTPCcalibDButil::ConnectCalPadTrees(baseDir,pattern,0); //create tree and attach calibration as friends
3335 :
3336 : //
3337 : // === add the calibPulser trees ======================================
3338 0 : TString inputTreesCalPad = gSystem->GetFromPipe(Form("ls %s/%s",baseDir.Data(), pattern.Data()));
3339 0 : TObjArray *arrInputTreesCalPad = inputTreesCalPad.Tokenize("\n");
3340 : //
3341 0 : for (Int_t itree=0; itree<arrInputTreesCalPad->GetEntriesFast(); ++itree) {
3342 0 : TFile *fin2 = TFile::Open(arrInputTreesCalPad->At(itree)->GetName());
3343 0 : TTree *tin = (TTree*)fin2->Get("calPads");
3344 0 : gROOT->cd();
3345 0 : TString friendName=arrInputTreesCalPad->At(itree)->GetName();
3346 0 : friendName.ReplaceAll("//","/");
3347 0 : friendName.ReplaceAll(baseDir.Data(),"");
3348 0 : friendName.ReplaceAll("^/","");
3349 0 : friendName.ReplaceAll("/",".");
3350 0 : friendName.ReplaceAll(".tree.",".");
3351 0 : friendName.ReplaceAll(".root","");
3352 0 : printf("%s\n", friendName.Data());
3353 0 : ::Info("AliTPCcalibDButil::ConnectCalPadTrees","%s",friendName.Data());
3354 0 : if (tMain==0) tMain=tin;
3355 0 : tMain->AddFriend(tin,friendName.Data());
3356 0 : TObjArray * branches=tin->GetListOfBranches();
3357 0 : Int_t nBranches=branches->GetEntries();
3358 0 : for (Int_t ibranch=0; ibranch<nBranches; ibranch++){
3359 0 : TString bname=branches->At(ibranch)->GetName();
3360 0 : if (bname.Contains(".")>0){
3361 0 : bname.ReplaceAll(".","");
3362 : // replace elements
3363 0 : tin->SetAlias((bname).Data(), (bname+".fElements").Data());
3364 0 : tMain->SetAlias((friendName+"."+bname).Data(), (friendName+"."+bname+".fElements").Data());
3365 : //
3366 : // make normalized values per chamber
3367 : //
3368 0 : if (branches->FindObject(bname+"_LTM")!=0){
3369 0 : tMain->SetAlias((friendName+"."+bname+"_MeanRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_Mean").Data());
3370 0 : tMain->SetAlias((friendName+"."+bname+"_MedianRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_Median").Data());
3371 0 : tMain->SetAlias((friendName+"."+bname+"_LTMRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_LTM").Data());
3372 0 : tMain->SetAlias((friendName+"."+bname+"_MeanDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_Mean").Data());
3373 0 : tMain->SetAlias((friendName+"."+bname+"_MedianDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_Median").Data());
3374 0 : tMain->SetAlias((friendName+"."+bname+"_LTMDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_LTM").Data());
3375 0 : }
3376 : // if (branches->FindObject(bname+"_Med3.")!=0){
3377 0 : tMain->SetAlias((friendName+"."+bname+"_Med3Ratio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Med3.fElements").Data());
3378 0 : tMain->SetAlias((friendName+"."+bname+"_Med5Ratio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Med5.fElements").Data());
3379 0 : tMain->SetAlias((friendName+"."+bname+"_Par6GRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Par6G.fElements").Data());
3380 0 : tMain->SetAlias((friendName+"."+bname+"_Med3Delta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Med3.fElements").Data());
3381 0 : tMain->SetAlias((friendName+"."+bname+"_Med5Delta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Med5.fElements").Data());
3382 0 : tMain->SetAlias((friendName+"."+bname+"_Par6GDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Par6G.fElements").Data());
3383 : // }
3384 0 : }
3385 0 : }
3386 0 : }
3387 : //
3388 : //
3389 : //
3390 : if (checkAliases){
3391 : // to be implemented
3392 : }
3393 : return tMain;
3394 0 : }
3395 :
3396 :
3397 :
|
