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