Constant extrapolation for gain calibration
[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
480//_____________________________________________________________________________________
abb20887 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
564//_____________________________________________________________________________________
892226be 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()
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)
875{
876 //
892226be 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 {
668e6ebb 1231 fRefDataQA=(AliTPCdataQA*)fRefDataQA->Clone();
949d8707 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;
58affc42 1356 return (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;
58affc42 1366 return (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 //
78f17711 1913 // Get the correction of the drift velocity using the offline laser tracks calbration
a23ba1c3 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);
78f17711 1921
1922 return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1923}
1924
1925Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1926 //
1927 // Get the correction of the drift velocity using the online laser tracks calbration
1928 //
1929 // run - run number
1930 // timeStamp - tim stamp in seconds
1931 // deltaT - integration period to calculate time0 offset
1932 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1933 // Note in case no data form both A and C side - the value from active side used
1934 TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1935
1936 Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1937 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1938 if (!param) return 0;
1939
1940 //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1941 dv*=param->GetDriftV()/2.61301900000000000e+06;
1942 if (dv>1e-20) dv=1/dv-1;
1943 else return 0;
1944 // T/P correction
1945 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1946
1947 AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
2cd3a697 1948 AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
78f17711 1949
1950 Double_t corrPTA=0;
1951 Double_t corrPTC=0;
1952
1953 if (temp&&press) {
1954 AliTPCCalibVdrift corr(temp,press,0);
1955 corrPTA=corr.GetPTRelative(timeStamp,0);
1956 corrPTC=corr.GetPTRelative(timeStamp,1);
1957 }
1958
1959 if (side==0) dv -= corrPTA;
1960 if (side==1) dv -= corrPTC;
1961 if (side==2) dv -= (corrPTA+corrPTC)/2;
1962
1963 return dv;
1964}
1965
1966Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1967 Int_t side, TObjArray * const array){
1968 //
1969 // common drift velocity retrieval for online and offline method
1970 //
a23ba1c3 1971 TGraphErrors *grlaserA=0;
1972 TGraphErrors *grlaserC=0;
1973 Double_t vlaserA=0, vlaserC=0;
1974 if (!array) return 0;
1975 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1976 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1977 Double_t deltaY;
5647625c 1978 if (grlaserA && grlaserA->GetN()>0) {
a23ba1c3 1979 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1980 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1981 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1982 }
5647625c 1983 if (grlaserC && grlaserC->GetN()>0) {
a23ba1c3 1984 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1985 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1986 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1987 }
1988 if (side==0) return vlaserA;
1989 if (side==1) return vlaserC;
1990 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1991 if (!grlaserA) return vlaserC;
1992 if (!grlaserC) return vlaserA;
1993 return mdrift;
1994}
1995
1996
1997Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1998 //
1999 // Get the correction of the drift velocity using the CE laser data
2000 // combining information from the CE, laser track calibration
2001 // and P/T calibration
2002 //
2003 // run - run number
2004 // timeStamp - tim stamp in seconds
2005 // deltaT - integration period to calculate time0 offset
2006 // side - 0 - A side, 1 - C side, 2 - mean from both sides
2007 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
0d1b4cf8 2008 if (!arrT) return 0;
a23ba1c3 2009 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
2010 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2011 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2012 //
2013 //
2014 Double_t corrPTA = 0, corrPTC=0;
2015 Double_t ltime0A = 0, ltime0C=0;
2016 Double_t gry=0;
2017 Double_t corrA=0, corrC=0;
2018 Double_t timeA=0, timeC=0;
5647625c 2019 const Double_t kEpsilon = 0.00001;
a23ba1c3 2020 TGraph *graphA = (TGraph*)arrT->At(72);
2021 TGraph *graphC = (TGraph*)arrT->At(73);
2022 if (!graphA && !graphC) return 0.;
2023 if (graphA &&graphA->GetN()>0) {
2024 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2025 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2026 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2027 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2028 if(ltime0A < kEpsilon) return 0;
a23ba1c3 2029 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2030 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2031 corrA-=corrPTA;
2032 }
2033 if (graphC&&graphC->GetN()>0){
2034 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2035 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2036 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2037 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2038 if(ltime0C < kEpsilon) return 0;
2039if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2040 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2041 corrC-=corrPTC;
2042 }
2043
2044 if (side ==0 ) return corrA;
2045 if (side ==1 ) return corrC;
2046 Double_t corrM= (corrA+corrC)*0.5;
2047 if (!graphA) corrM=corrC;
2048 if (!graphC) corrM=corrA;
2049 return corrM;
2050}
2051
cc65e4f5 2052Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2053 //
2054 // return drift velocity using the TPC-ITS matchin method
2055 // return also distance to the closest point
2056 //
2057 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2058 TGraphErrors *graph=0;
2059 dist=0;
2060 if (!array) return 0;
5647625c 2061 //array->ls();
cc65e4f5 2062 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
2063 if (!graph) return 0;
2064 Double_t deltaY;
2065 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2066 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2067 return value;
2068}
2069
2070Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2071 //
2072 // Get time dependent time 0 (trigger delay in cm) correction
2073 // Arguments:
2074 // timestamp - timestamp
2075 // run - run number
2076 //
2077 // Notice - Extrapolation outside of calibration range - using constant function
2078 //
2079 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2080 TGraphErrors *graph=0;
2081 dist=0;
2082 if (!array) return 0;
2083 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
2084 if (!graph) return 0;
2085 Double_t deltaY;
2086 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2087 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2088 return value;
2089}
2090
2091
a23ba1c3 2092
2093
2094
2095Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2096 //
2097 // VERY obscure method - we need something in framework
2098 // Find the TPC runs with temperature OCDB entry
2099 // cache the start and end of the run
2100 //
2101 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2102 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2103 if (!storage) return 0;
2104 TString path=storage->GetURI();
2105 TString runsT;
2106 {
2107 TString command;
2108 if (path.Contains("local")){ // find the list if local system
2109 path.ReplaceAll("local://","");
2110 path+="TPC/Calib/Temperature";
2111 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2112 }
2113 runsT=gSystem->GetFromPipe(command);
2114 }
2115 TObjArray *arr= runsT.Tokenize("r");
2116 if (!arr) return 0;
2117 //
2118 TArrayI indexes(arr->GetEntries());
2119 TArrayI runs(arr->GetEntries());
2120 Int_t naccept=0;
2121 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2122 Int_t irunN = atoi(arr->At(irun)->GetName());
2123 if (irunN<startRun) continue;
2124 if (irunN>stopRun) continue;
2125 runs[naccept]=irunN;
2126 naccept++;
2127 }}
2128 fRuns.Set(naccept);
2129 fRunsStart.Set(fRuns.fN);
2130 fRunsStop.Set(fRuns.fN);
2131 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2132 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2133
2134 //
2135 AliCDBEntry * entry = 0;
2136 {for (Int_t irun=0;irun<fRuns.fN; irun++){
2137 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2138 if (!entry) continue;
2139 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2140 if (!tmpRun) continue;
2141 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2142 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2143 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
2144 }}
2145 return fRuns.fN;
2146}
2147
2148
2149Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2150 //
2151 // binary search - find the run for given time stamp
2152 //
2153 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2154 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2155 Int_t cindex = -1;
2156 for (Int_t index=index0; index<=index1; index++){
2157 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2158 if (debug) {
2159 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2160 }
2161 }
2162 if (cindex<0) cindex =(index0+index1)/2;
2163 if (cindex<0) {
2164 return 0;
2165 }
2166 return fRuns[cindex];
2167}
2168
2169
2170
2171
2172
2173TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2174 //
2175 // filter outlyer measurement
2176 // Only points around median +- sigmaCut filtered
2177 if (!graph) return 0;
2178 Int_t kMinPoints=2;
2179 Int_t npoints0 = graph->GetN();
2180 Int_t npoints=0;
2181 Float_t rmsY=0;
a23ba1c3 2182 //
2183 //
2184 if (npoints0<kMinPoints) return 0;
a9f487a1 2185
2186 Double_t *outx=new Double_t[npoints0];
2187 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2188 for (Int_t iter=0; iter<3; iter++){
2189 npoints=0;
2190 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2191 if (graph->GetY()[ipoint]==0) continue;
2192 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2193 outx[npoints] = graph->GetX()[ipoint];
2194 outy[npoints] = graph->GetY()[ipoint];
2195 npoints++;
2196 }
2197 if (npoints<=1) break;
2198 medianY =TMath::Median(npoints,outy);
2199 rmsY =TMath::RMS(npoints,outy);
2200 }
2201 TGraph *graphOut=0;
2202 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2203 delete [] outx;
2204 delete [] outy;
a23ba1c3 2205 return graphOut;
2206}
2207
2208
2209TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2210 //
2211 // filter outlyer measurement
2212 // Only points around median +- cut filtered
2213 if (!graph) return 0;
2214 Int_t kMinPoints=2;
2215 Int_t npoints0 = graph->GetN();
2216 Int_t npoints=0;
2217 Float_t rmsY=0;
a23ba1c3 2218 //
2219 //
2220 if (npoints0<kMinPoints) return 0;
a9f487a1 2221
2222 Double_t *outx=new Double_t[npoints0];
2223 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2224 for (Int_t iter=0; iter<3; iter++){
2225 npoints=0;
2226 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2227 if (graph->GetY()[ipoint]==0) continue;
2228 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2229 outx[npoints] = graph->GetX()[ipoint];
2230 outy[npoints] = graph->GetY()[ipoint];
2231 npoints++;
2232 }
2233 if (npoints<=1) break;
2234 medianY =TMath::Median(npoints,outy);
2235 rmsY =TMath::RMS(npoints,outy);
2236 }
2237 TGraph *graphOut=0;
2238 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2239 delete [] outx;
2240 delete [] outy;
a23ba1c3 2241 return graphOut;
2242}
2243
2244
2245
8166d5d7 2246TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
a23ba1c3 2247 //
2248 // filter outlyer measurement
2249 // Only points with normalized errors median +- sigmaCut filtered
2250 //
2251 Int_t kMinPoints=10;
2252 Int_t npoints0 = graph->GetN();
2253 Int_t npoints=0;
2254 Float_t medianErr=0, rmsErr=0;
a9f487a1 2255 //
2256 //
2257 if (npoints0<kMinPoints) return 0;
2258
a23ba1c3 2259 Double_t *outx=new Double_t[npoints0];
2260 Double_t *outy=new Double_t[npoints0];
2261 Double_t *erry=new Double_t[npoints0];
2262 Double_t *nerry=new Double_t[npoints0];
2263 Double_t *errx=new Double_t[npoints0];
a9f487a1 2264
a23ba1c3 2265 for (Int_t iter=0; iter<3; iter++){
2266 npoints=0;
2267 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2268 nerry[npoints] = graph->GetErrorY(ipoint);
2269 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2270 erry[npoints] = graph->GetErrorY(ipoint);
2271 outx[npoints] = graph->GetX()[ipoint];
2272 outy[npoints] = graph->GetY()[ipoint];
2273 errx[npoints] = graph->GetErrorY(ipoint);
2274 npoints++;
2275 }
2276 if (npoints==0) break;
2277 medianErr=TMath::Median(npoints,erry);
2278 medianY =TMath::Median(npoints,outy);
2279 rmsErr =TMath::RMS(npoints,erry);
2280 }
2281 TGraphErrors *graphOut=0;
2282 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2283 delete []outx;
2284 delete []outy;
a23ba1c3 2285 delete []erry;
a9f487a1 2286 delete []nerry;
2287 delete []errx;
a23ba1c3 2288 return graphOut;
2289}
2290
2291
2292void AliTPCcalibDButil::Sort(TGraph *graph){
2293 //
2294 // sort array - neccessay for approx
2295 //
2296 Int_t npoints = graph->GetN();
2297 Int_t *indexes=new Int_t[npoints];
2298 Double_t *outx=new Double_t[npoints];
2299 Double_t *outy=new Double_t[npoints];
2300 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2301 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2302 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2303 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2304 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2305
a9f487a1 2306 delete [] indexes;
2307 delete [] outx;
2308 delete [] outy;
a23ba1c3 2309}
2310void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2311 //
2312 // smmoth graph - mean on the interval
2313 //
2314 Sort(graph);
2315 Int_t npoints = graph->GetN();
2316 Double_t *outy=new Double_t[npoints];
0d1b4cf8 2317
a23ba1c3 2318 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2319 Double_t lx=graph->GetX()[ipoint];
2320 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2321 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2322 if (index0<0) index0=0;
2323 if (index1>=npoints-1) index1=npoints-1;
2324 if ((index1-index0)>1){
2325 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2326 }else{
2327 outy[ipoint]=graph->GetY()[ipoint];
2328 }
2329 }
0d1b4cf8 2330 // TLinearFitter fitter(3,"pol2");
2331// for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2332// Double_t lx=graph->GetX()[ipoint];
2333// Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2334// Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2335// if (index0<0) index0=0;
2336// if (index1>=npoints-1) index1=npoints-1;
2337// fitter.ClearPoints();
2338// for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2339// if ((index1-index0)>1){
2340// outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2341// }else{
2342// outy[ipoint]=graph->GetY()[ipoint];
2343// }
2344// }
2345
2346
2347
a23ba1c3 2348 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2349 graph->GetY()[ipoint] = outy[ipoint];
2350 }
2351 delete[] outy;
2352}
2353
8166d5d7 2354Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
a23ba1c3 2355 //
2356 // Use constant interpolation outside of range
2357 //
2358 if (!graph) {
2359 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2360 return 0;
2361 }
5647625c 2362
a23ba1c3 2363 if (graph->GetN()<1){
5647625c 2364 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
a23ba1c3 2365 return 0;
2366 }
5647625c 2367
2368
a23ba1c3 2369 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2370 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
5647625c 2371
2b68ab9c 2372 // printf("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref));
5647625c 2373
2374 if(graph->GetN()==1)
2375 return graph->Eval(graph->GetX()[0]);
2376
2377
2378 return graph->Eval(xref);
a23ba1c3 2379}
2380
108a67db 2381Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2382 //
2383 // Use constant interpolation outside of range also for spline fits
2384 //
2385 if (!graph) {
2386 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2387 return 0;
2388 }
2389 if (graph->GetKnots()<1){
2390 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2391 return 0;
2392 }
2393 if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2394 if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2395 return graph->Eval( xref);
2396}
2397
1e722a63 2398Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2399 //
2400 // Filter DCS sensor information
2401 // ymin - minimal value
2402 // ymax - max value
2403 // maxdy - maximal deirivative
2404 // sigmaCut - cut on values and derivative in terms of RMS distribution
2405 // Return value - accepted fraction
2406 //
2407 // Algorithm:
2408 //
2409 // 0. Calculate median and rms of values in specified range
2410 // 1. Filter out outliers - median+-sigmaCut*rms
2411 // values replaced by median
2412 //
2413 AliSplineFit * fit = sensor->GetFit();
2414 if (!fit) return 0.;
2415 Int_t nknots = fit->GetKnots();
2416 if (nknots==0) {
2417 delete fit;
2418 sensor->SetFit(0);
2419 return 0;
2420 }
2421 //
2422 Double_t *yin0 = new Double_t[nknots];
2423 Double_t *yin1 = new Double_t[nknots];
2424 Int_t naccept=0;
2425
2426 for (Int_t iknot=0; iknot< nknots; iknot++){
2427 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2428 yin0[naccept] = fit->GetY0()[iknot];
2429 yin1[naccept] = fit->GetY1()[iknot];
2430 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2431 naccept++;
2432 }
2433 }
2434 if (naccept<1) {
2435 delete fit;
2436 sensor->SetFit(0);
a9f487a1 2437 delete [] yin0;
2438 delete [] yin1;
1e722a63 2439 return 0.;
2440 }
0d1b4cf8 2441
1e722a63 2442 Double_t medianY0=0, medianY1=0;
2443 Double_t rmsY0 =0, rmsY1=0;
2444 medianY0 = TMath::Median(naccept, yin0);
2445 medianY1 = TMath::Median(naccept, yin1);
2446 rmsY0 = TMath::RMS(naccept, yin0);
2447 rmsY1 = TMath::RMS(naccept, yin1);
2448 naccept=0;
2449 //
2450 // 1. Filter out outliers - median+-sigmaCut*rms
2451 // values replaced by median
2452 // if replaced the derivative set to 0
2453 //
2454 for (Int_t iknot=0; iknot< nknots; iknot++){
2455 Bool_t isOK=kTRUE;
2456 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2457 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2458 if (nknots<2) fit->GetY1()[iknot]=0;
2459 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2460 if (!isOK){
2461 fit->GetY0()[iknot]=medianY0;
2462 fit->GetY1()[iknot]=0;
2463 }else{
2464 naccept++;
2465 }
2466 }
2467 delete [] yin0;
2468 delete [] yin1;
2469 return Float_t(naccept)/Float_t(nknots);
2470}
2471
2472Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2473 //
2474 // Filter temperature array
2475 // tempArray - array of temperatures -
2476 // ymin - minimal accepted temperature - default 15
2477 // ymax - maximal accepted temperature - default 22
2478 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2479 // return value - fraction of filtered sensors
2480 const Double_t kMaxDy=0.1;
2481 Int_t nsensors=tempArray->NumSensors();
2482 if (nsensors==0) return 0.;
2483 Int_t naccept=0;
2484 for (Int_t isensor=0; isensor<nsensors; isensor++){
2485 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2486 if (!sensor) continue;
2487 //printf("%d\n",isensor);
2488 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2489 if (sensor->GetFit()==0){
0d1b4cf8 2490 //delete sensor;
1e722a63 2491 tempArray->RemoveSensorNum(isensor);
2492 }else{
2493 naccept++;
2494 }
2495 }
2496 return Float_t(naccept)/Float_t(nsensors);
2497}
2498
a23ba1c3 2499
8166d5d7 2500void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2501 //
2502 // Filter CE data
1e722a63 2503 // Input parameters:
2504 // deltaT - smoothing window (in seconds)
2505 // cutAbs - max distance of the time info to the median (in time bins)
2506 // cutSigma - max distance (in the RMS)
2507 // pcstream - optional debug streamer to store original and filtered info
2508 // Hardwired parameters:
2509 // kMinPoints =10; // minimal number of points to define the CE
2510 // kMinSectors=12; // minimal number of sectors to define sideCE
2511 // Algorithm:
2512 // 0. Filter almost emty graphs (kMinPoints=10)
2513 // 1. calculate median and RMS per side
2514 // 2. Filter graphs - in respect with side medians
2515 // - cutAbs and cutDelta used
2516 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2517 // 4. Calculate mean for A side and C side
2518 //
2519 const Int_t kMinPoints =10; // minimal number of points to define the CE
2520 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2521 const Int_t kMinTime =400; // minimal arrival time of CE
2522 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
a23ba1c3 2523 Double_t medianY=0;
1e722a63 2524 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
a23ba1c3 2525 if (!cearray) return;
1e722a63 2526 Double_t tmin=-1;
2527 Double_t tmax=-1;
2528 //
2529 //
a23ba1c3 2530 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2cd3a697 2531 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
a23ba1c3 2532 if ( tempMapCE && cavernPressureCE){
1e722a63 2533 //
c5dfa2dc 2534 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2535 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2536 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2537 Bool_t isOK=kTRUE;
1e722a63 2538 if (isOK) {
2539 // recalculate P/T correction map for time of the CE
2540 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2541 driftCalib->SetName("driftPTCE");
2542 driftCalib->SetTitle("driftPTCE");
2543 cearray->AddLast(driftCalib);
2544 }
a23ba1c3 2545 }
1e722a63 2546 //
2547 // 0. Filter almost emty graphs
2548 //
a23ba1c3 2549
1e722a63 2550 for (Int_t i=0; i<72;i++){
a23ba1c3 2551 TGraph *graph= (TGraph*)arrT->At(i);
29797f44 2552 if (!graph) continue;
2553 graph->Sort();
a23ba1c3 2554 if (graph->GetN()<kMinPoints){
2555 arrT->AddAt(0,i);
2556 delete graph; // delete empty graph
2557 continue;
2558 }
1e722a63 2559 if (tmin<0) tmin = graph->GetX()[0];
2560 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2561 //
2562 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2563 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2564 }
2565 //
2566 // 1. calculate median and RMS per side
2567 //
2568 TArrayF arrA(100000), arrC(100000);
2569 Int_t nA=0, nC=0;
2570 Double_t medianA=0, medianC=0;
2571 Double_t rmsA=0, rmsC=0;
2572 for (Int_t isec=0; isec<72;isec++){
2573 TGraph *graph= (TGraph*)arrT->At(isec);
2574 if (!graph) continue;
2575 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2576 if (graph->GetY()[ipoint]<kMinTime) continue;
2577 if (nA>=arrA.fN) arrA.Set(nA*2);
2578 if (nC>=arrC.fN) arrC.Set(nC*2);
2579 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2580 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2581 }
2582 }
2583 if (nA>0){
2584 medianA=TMath::Median(nA,arrA.fArray);
2585 rmsA =TMath::RMS(nA,arrA.fArray);
2586 }
2587 if (nC>0){
2588 medianC=TMath::Median(nC,arrC.fArray);
2589 rmsC =TMath::RMS(nC,arrC.fArray);
2590 }
2591 //
2592 // 2. Filter graphs - in respect with side medians
2593 //
2594 TArrayD vecX(100000), vecY(100000);
2595 for (Int_t isec=0; isec<72;isec++){
2596 TGraph *graph= (TGraph*)arrT->At(isec);
2597 if (!graph) continue;
2598 Double_t median = (isec%36<18) ? medianA: medianC;
2599 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2600 Int_t naccept=0;
29797f44 2601 // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2602 for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
1e722a63 2603 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2604 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2605 vecX[naccept]= graph->GetX()[ipoint];
2606 vecY[naccept]= graph->GetY()[ipoint];
2607 naccept++;
2608 }
2609 if (naccept<kMinPoints){
2610 arrT->AddAt(0,isec);
2611 delete graph; // delete empty graph
2612 continue;
2613 }
2614 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2615 delete graph;
2616 arrT->AddAt(graph2,isec);
2617 }
2618 //
2619 // 3. Cut in respect wit the graph median
2620 //
2621 for (Int_t i=0; i<72;i++){
2622 TGraph *graph= (TGraph*)arrT->At(i);
2623 if (!graph) continue;
2624 //
2625 // filter in range
2626 //
2627 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
a23ba1c3 2628 if (!graphTS0) continue;
2629 if (graphTS0->GetN()<kMinPoints) {
2630 delete graphTS0;
2631 delete graph;
2632 arrT->AddAt(0,i);
2633 continue;
2634 }
1e722a63 2635 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
a23ba1c3 2636 graphTS->Sort();
2637 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2638 if (pcstream){
2639 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2640 (*pcstream)<<"filterCE"<<
2641 "run="<<run<<
2642 "isec="<<i<<
2643 "mY="<<medianY<<
2644 "graph.="<<graph<<
2645 "graphTS0.="<<graphTS0<<
2646 "graphTS.="<<graphTS<<
2647 "\n";
2648 }
2649 delete graphTS0;
2650 if (!graphTS) continue;
2651 arrT->AddAt(graphTS,i);
2652 delete graph;
2653 }
1e722a63 2654 //
2655 // Recalculate the mean time A side C side
2656 //
2657 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2658 Int_t meanPoints=(nA+nC)/72; // mean number of points
2659 for (Int_t itime=0; itime<200; itime++){
2660 nA=0, nC=0;
2661 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2662 for (Int_t i=0; i<72;i++){
2663 TGraph *graph= (TGraph*)arrT->At(i);
2664 if (!graph) continue;
2665 if (graph->GetN()<(meanPoints/4)) continue;
2666 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2667 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2668 }
2669 xA[itime]=time;
2670 xC[itime]=time;
2671 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2672 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2673 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2674 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2675 }
2676 //
2677 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2678 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2679 if (pcstream){
2680 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2681 (*pcstream)<<"filterAC"<<
2682 "run="<<run<<
2683 "nA="<<nA<<
2684 "nC="<<nC<<
2685 "rmsTA="<<rmsTA<<
2686 "rmsTC="<<rmsTC<<
2687 "\n";
2688 }
2689 //
2690 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2691 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2692 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2693 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2694 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2695 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2696 delete grA;
2697 delete grC;
2698 if (nA<kMinSectors) arrT->AddAt(0,72);
2699 else arrT->AddAt(graphTSA,72);
2700 if (nC<kMinSectors) arrT->AddAt(0,73);
2701 else arrT->AddAt(graphTSC,73);
a23ba1c3 2702}
2703
2704
8166d5d7 2705void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2706 //
2707 // Filter Drift velocity measurement using the tracks
2708 // 0. remove outlyers - error based
2709 // cutSigma
2710 //
2711 //
2712 const Int_t kMinPoints=1; // minimal number of points to define value
2713 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2714 Double_t medianY=0;
2715 if (!arrT) return;
2716 for (Int_t i=0; i<arrT->GetEntries();i++){
14301155 2717 TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
a23ba1c3 2718 if (!graph) continue;
2719 if (graph->GetN()<kMinPoints){
2720 delete graph;
2721 arrT->AddAt(0,i);
2722 continue;
2723 }
5647625c 2724 TGraphErrors *graph2 = NULL;
2725 if(graph->GetN()<10) {
2726 graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2727 if (!graph2) {
2728 delete graph; arrT->AddAt(0,i); continue;
2729 }
2730 }
2731 else {
2732 graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2733 if (!graph2) {
2734 delete graph; arrT->AddAt(0,i); continue;
2735 }
a23ba1c3 2736 }
2737 if (graph2->GetN()<1) {
2738 delete graph; arrT->AddAt(0,i); continue;
2739 }
2740 graph2->SetName(graph->GetName());
2741 graph2->SetTitle(graph->GetTitle());
2742 arrT->AddAt(graph2,i);
2743 if (pcstream){
2744 (*pcstream)<<"filterTracks"<<
2745 "run="<<run<<
2746 "isec="<<i<<
2747 "mY="<<medianY<<
2748 "graph.="<<graph<<
2749 "graph2.="<<graph2<<
2750 "\n";
2751 }
2752 delete graph;
2753 }
2754}
2755
2756
a23ba1c3 2757
2758
2759
2760Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2761 //
2762 //
2763 // get laser time offset
2764 // median around timeStamp+-deltaT
2765 // QA - chi2 needed for later usage - to be added
2766 // - currently cut on error
2767 //
2768 Int_t kMinPoints=1;
2769 Double_t kMinDelay=0.01;
2770 Double_t kMinDelayErr=0.0001;
2771
2772 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2773 if (!array) return 0;
2774 TGraphErrors *tlaser=0;
2775 if (array){
2776 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2777 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2778 }
2779 if (!tlaser) return 0;
2780 Int_t npoints0= tlaser->GetN();
2781 if (npoints0==0) return 0;
2782 Double_t *xlaser = new Double_t[npoints0];
2783 Double_t *ylaser = new Double_t[npoints0];
2784 Int_t npoints=0;
2785 for (Int_t i=0;i<npoints0;i++){
2786 //printf("%d\n",i);
2787 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2788 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2789 xlaser[npoints]=tlaser->GetX()[npoints];
2790 ylaser[npoints]=tlaser->GetY()[npoints];
2791 npoints++;
2792 }
2793 //
2794 //
2795 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2796 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2797 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2798 if (index0<0) index0=0;
2799 if (index1>=npoints-1) index1=npoints-1;
2800 if (index1-index0<kMinPoints) return 0;
2801 //
2802 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2803 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2804 delete [] ylaser;
2805 delete [] xlaser;
2806 return mean;
2807}
0d1b4cf8 2808
2809
2810
2811
8166d5d7 2812void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
0d1b4cf8 2813 //
2814 // Filter Goofie data
1fabc823 2815 // goofieArray - points will be filtered
2816 // deltaT - smmothing time window
2817 // cutSigma - outler sigma cut in rms
2818 // minVn, maxVd- range absolute cut for variable vd/pt
2819 // - to be tuned
0d1b4cf8 2820 //
0d1b4cf8 2821 // Ignore goofie if not enough points
2822 //
2823 const Int_t kMinPoints = 3;
2824 //
2825
2826 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2827 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2828 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2829 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2830 if (!graphvd) return;
2831 if (graphvd->GetN()<kMinPoints){
2832 delete graphvd;
2833 goofieArray->GetSensorNum(2)->SetGraph(0);
2834 return;
2835 }
2836 //
2837 // 1. Caluclate medians of critical variables
2838 // drift velcocity
2839 // P/T
2840 // area near peak
2841 // area far peak
2842 //
2843 Double_t medianpt=0;
1fabc823 2844 Double_t medianvd=0, sigmavd=0;
0d1b4cf8 2845 Double_t medianan=0;
2846 Double_t medianaf=0;
1fabc823 2847 Int_t entries=graphvd->GetN();
2848 Double_t yvdn[10000];
2849 Int_t nvd=0;
2850 //
2851 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2852 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2853 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2854 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2855 yvdn[nvd++]=graphvd->GetY()[ipoint];
2856 }
2857 if (nvd<kMinPoints){
2858 delete graphvd;
2859 goofieArray->GetSensorNum(2)->SetGraph(0);
2860 return;
2861 }
2862 //
2863 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2864 if (nuni>=kMinPoints){
2865 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2866 }else{
2867 medianvd = TMath::Median(nvd, yvdn);
2868 }
2869
0d1b4cf8 2870 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2871 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2872 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2873 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2874 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2875 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
0d1b4cf8 2876 delete graphpt0;
2877 delete graphpt1;
2878 delete graphan0;
2879 delete graphan1;
2880 delete graphaf0;
2881 delete graphaf1;
2882 //
2883 // 2. Make outlyer graph
2884 //
1fabc823 2885 Int_t nOK=0;
0d1b4cf8 2886 TGraph graphOut(*graphvd);
2887 for (Int_t i=0; i<entries;i++){
2888 //
2889 Bool_t isOut=kFALSE;
1fabc823 2890 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2891 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2892
2893 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2894 isOut|=kTRUE;
0d1b4cf8 2895 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
1fabc823 2896 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2897 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
0d1b4cf8 2898 graphOut.GetY()[i]= (isOut)?1:0;
1fabc823 2899 if (!isOut) nOK++;
0d1b4cf8 2900 }
1fabc823 2901 if (nOK<kMinPoints) {
2902 delete graphvd;
2903 goofieArray->GetSensorNum(2)->SetGraph(0);
2904 return;
2905 }
0d1b4cf8 2906 //
2907 // 3. Filter out outlyers - and smooth
2908 //
2909 TVectorF vmedianArray(goofieArray->NumSensors());
2910 TVectorF vrmsArray(goofieArray->NumSensors());
2911 Double_t xnew[10000];
2912 Double_t ynew[10000];
2913 TObjArray junk;
2914 junk.SetOwner(kTRUE);
2915 Bool_t isOK=kTRUE;
2916 //
2917 //
2918 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2919 isOK=kTRUE;
2920 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2921 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2922 //
2923 if (!sensor) continue;
2924 graphOld = sensor->GetGraph();
2925 if (graphOld) {
2926 sensor->SetGraph(0);
2927 Int_t nused=0;
2928 for (Int_t i=0;i<entries;i++){
2929 if (graphOut.GetY()[i]>0.5) continue;
2930 xnew[nused]=graphOld->GetX()[i];
2931 ynew[nused]=graphOld->GetY()[i];
2932 nused++;
2933 }
2934 graphNew = new TGraph(nused,xnew,ynew);
2935 junk.AddLast(graphNew);
2936 junk.AddLast(graphOld);
2937 Double_t median=0;
2938 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2939 if (graphNew0!=0){
2940 junk.AddLast(graphNew0);
2941 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2942 if (graphNew1!=0){
1fabc823 2943 junk.AddLast(graphNew1);
0d1b4cf8 2944 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2945 if (graphNew2!=0) {
1fabc823 2946 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
0d1b4cf8 2947 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2948 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2949 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
1fabc823 2950 printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
0d1b4cf8 2951 vmedianArray[isensor]=median;
0d1b4cf8 2952 //
2953 }
2954 }
2955 }
2956 }
2957 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2958 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2959 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2960 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2961 (*pcstream)<<"goofieA"<<
2962 Form("isOK_%d.=",isensor)<<isOK<<
2963 Form("s_%d.=",isensor)<<sensor<<
2964 Form("gr_%d.=",isensor)<<graphOld<<
2965 Form("gr0_%d.=",isensor)<<graphNew0<<
2966 Form("gr1_%d.=",isensor)<<graphNew1<<
2967 Form("gr2_%d.=",isensor)<<graphNew2;
1fabc823 2968 if (isOK) sensor->SetGraph(graphNew2);
2969 }
2970 (*pcstream)<<"goofieA"<<
2971 "vmed.="<<&vmedianArray<<
2972 "vrms.="<<&vrmsArray<<
2973 "\n";
2974 junk.Delete(); // delete temoprary graphs
0d1b4cf8 2975
2976}
2977
2978
a980538f 2979
2980
2981
8166d5d7 2982TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
a980538f 2983 //
2984 // Make a statistic matrix
2985 // Input parameters:
2986 // array - TObjArray of AliRelKalmanAlign
2987 // minFraction - minimal ration of accepted tracks
2988 // minStat - minimal statistic (number of accepted tracks)
2989 // maxvd - maximal deviation for the 1
2990 // Output matrix:
2991 // columns - Mean, Median, RMS
2992 // row - parameter type (rotation[3], translation[3], drift[3])
2993 if (!array) return 0;
2994 if (array->GetEntries()<=0) return 0;
cc65e4f5 2995 // Int_t entries = array->GetEntries();
a980538f 2996 Int_t entriesFast = array->GetEntriesFast();
2997 TVectorD state(9);
2998 TVectorD *valArray[9];
2999 for (Int_t i=0; i<9; i++){
3000 valArray[i] = new TVectorD(entriesFast);
3001 }
3002 Int_t naccept=0;
3003 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3004 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3005 if (!kalman) continue;
3006 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3007 if (kalman->GetNUpdates()<minStat) continue;
3008 if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
3009 kalman->GetState(state);
3010 for (Int_t ipar=0; ipar<9; ipar++)
3011 (*valArray[ipar])[naccept]=state[ipar];
3012 naccept++;
3013 }
5647625c 3014 //if (naccept<2) return 0;
3015 if (naccept<1) return 0;
a980538f 3016 TMatrixD *pstat=new TMatrixD(9,3);
3017 TMatrixD &stat=*pstat;
3018 for (Int_t ipar=0; ipar<9; ipar++){
3019 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3020 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3021 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3022 }
3023 return pstat;
3024}
3025
3026
8166d5d7 3027TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
a980538f 3028 //
3029 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3030 // Input:
3031 // array - input array
3032 // stat - mean parameters statistic
3033 // direction -
3034 // sigmaCut - maximal allowed deviation from mean in terms of RMS
3035 if (!array) return 0;
3036 if (array->GetEntries()<=0) return 0;
32100b1d 3037 if (!(&stat)) return 0;
cc65e4f5 3038 // error increase in 1 hour
3039 const Double_t kerrsTime[9]={
3040 0.00001, 0.00001, 0.00001,
3041 0.001, 0.001, 0.001,
dfb83639 3042 0.002, 0.01, 0.001};
cc65e4f5 3043 //
3044 //
a980538f 3045 Int_t entries = array->GetEntriesFast();
3046 TObjArray *sArray= new TObjArray(entries);
3047 AliRelAlignerKalman * sKalman =0;
3048 TVectorD state(9);
3049 for (Int_t i=0; i<entries; i++){
3050 Int_t index=(direction)? entries-i-1:i;
3051 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3052 if (!kalman) continue;
3053 Bool_t isOK=kTRUE;
3054 kalman->GetState(state);
3055 for (Int_t ipar=0; ipar<9; ipar++){
3056 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3057 }
3058 if (!sKalman &&isOK) {
3059 sKalman=new AliRelAlignerKalman(*kalman);
3060 sKalman->SetRejectOutliers(kFALSE);
3061 sKalman->SetRunNumber(kalman->GetRunNumber());
3062 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3063 }
3064 if (!sKalman) continue;
3065 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
cc65e4f5 3066 for (Int_t ipar=0; ipar<9; ipar++){
3067// (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3068// (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3069// (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3070 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3071 }
a980538f 3072 sKalman->SetRunNumber(kalman->GetRunNumber());
3073 if (!isOK) sKalman->SetRunNumber(0);
3074 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3075 if (!isOK) continue;
3076 sKalman->SetRejectOutliers(kFALSE);
3077 sKalman->SetRunNumber(kalman->GetRunNumber());
3078 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3079 sKalman->Merge(kalman);
3080 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3081 //sKalman->Print();
3082 }
3083 return sArray;
3084}
3085
8166d5d7 3086TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
cc65e4f5 3087 //
3088 // Merge 2 RelKalman arrays
3089 // Input:
3090 // arrayP - rel kalman in direction plus
3091 // arrayM - rel kalman in direction minus
3092 if (!arrayP) return 0;
3093 if (arrayP->GetEntries()<=0) return 0;
3094 if (!arrayM) return 0;
3095 if (arrayM->GetEntries()<=0) return 0;
3a4b0aed 3096
cc65e4f5 3097 Int_t entries = arrayP->GetEntriesFast();
3098 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3a4b0aed 3099
cc65e4f5 3100 for (Int_t i=0; i<entries; i++){
3101 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3102 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3103 if (!kalmanP) continue;
3104 if (!kalmanM) continue;
3a4b0aed 3105
3106 AliRelAlignerKalman *kalman = NULL;
3107 if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3108 kalman = new AliRelAlignerKalman(*kalmanP);
3109 kalman->Merge(kalmanM);
3110 }
3111 else if (kalmanP->GetRunNumber() == 0) {
3112 kalman = new AliRelAlignerKalman(*kalmanM);
3113 }
3114 else if (kalmanM->GetRunNumber() == 0) {
3115 kalman = new AliRelAlignerKalman(*kalmanP);
3116 }
3117 else
3118 continue;
3119
cc65e4f5 3120 array->AddAt(kalman,i);
3121 }
3a4b0aed 3122
cc65e4f5 3123 return array;
3124}