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