]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Base/AliTPCcalibDButil.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Extracting Branch and Revision from Git.
[u/mrichter/AliRoot.git] / TPC / Base / 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>
df5d3222 43#include <AliSplineFit.h>
892226be 44#include "AliTPCcalibDB.h"
45#include "AliTPCCalPad.h"
46#include "AliTPCCalROC.h"
47#include "AliTPCROC.h"
48#include "AliTPCmapper.h"
49#include "AliTPCParam.h"
6e7d7dc4 50#include "AliTPCCalibRaw.h"
949d8707 51#include "AliTPCPreprocessorOnline.h"
52#include "AliTPCdataQA.h"
a23ba1c3 53#include "AliLog.h"
892226be 54#include "AliTPCcalibDButil.h"
a23ba1c3 55#include "AliTPCCalibVdrift.h"
1fabc823 56#include "AliMathBase.h"
a980538f 57#include "AliRelAlignerKalman.h"
892226be 58
8166d5d7 59const Float_t kAlmost0=1.e-30;
60
892226be 61ClassImp(AliTPCcalibDButil)
62AliTPCcalibDButil::AliTPCcalibDButil() :
63 TObject(),
1e722a63 64 fCalibDB(0),
892226be 65 fPadNoise(0x0),
66 fPedestals(0x0),
67 fPulserTmean(0x0),
68 fPulserTrms(0x0),
69 fPulserQmean(0x0),
70 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
71 fCETmean(0x0),
72 fCETrms(0x0),
73 fCEQmean(0x0),
74 fALTROMasked(0x0),
6e7d7dc4 75 fCalibRaw(0x0),
949d8707 76 fDataQA(0x0),
77 fRefMap(0x0),
78 fCurrentRefMap(0x0),
79 fRefValidity(""),
7390f655 80 fRefPadNoise(0x0),
81 fRefPedestals(0x0),
949d8707 82 fRefPedestalMasked(0x0),
7390f655 83 fRefPulserTmean(0x0),
84 fRefPulserTrms(0x0),
85 fRefPulserQmean(0x0),
86 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
949d8707 87 fRefPulserMasked(0x0),
7390f655 88 fRefCETmean(0x0),
89 fRefCETrms(0x0),
90 fRefCEQmean(0x0),
949d8707 91 fRefCEMasked(0x0),
92 fRefALTROFPED(0x0),
93 fRefALTROZsThr(0x0),
94 fRefALTROAcqStart(0x0),
95 fRefALTROAcqStop(0x0),
7390f655 96 fRefALTROMasked(0x0),
97 fRefCalibRaw(0x0),
949d8707 98 fRefDataQA(0x0),
892226be 99 fGoofieArray(0x0),
100 fMapper(new AliTPCmapper(0x0)),
101 fNpulserOutliers(-1),
7390f655 102 fIrocTimeOffset(0),
892226be 103 fCETmaxLimitAbs(1.5),
104 fPulTmaxLimitAbs(1.5),
105 fPulQmaxLimitAbs(5),
a23ba1c3 106 fPulQminLimit(11),
107 fRuns(0), // run list with OCDB info
108 fRunsStart(0), // start time for given run
109 fRunsStop(0) // stop time for given run
892226be 110{
111 //
112 // Default ctor
113 //
114}
115//_____________________________________________________________________________________
116AliTPCcalibDButil::~AliTPCcalibDButil()
117{
118 //
119 // dtor
120 //
121 delete fPulserOutlier;
7390f655 122 delete fRefPulserOutlier;
892226be 123 delete fMapper;
7390f655 124 if (fRefPadNoise) delete fRefPadNoise;
125 if (fRefPedestals) delete fRefPedestals;
949d8707 126 if (fRefPedestalMasked) delete fRefPedestalMasked;
7390f655 127 if (fRefPulserTmean) delete fRefPulserTmean;
128 if (fRefPulserTrms) delete fRefPulserTrms;
129 if (fRefPulserQmean) delete fRefPulserQmean;
949d8707 130 if (fRefPulserMasked) delete fRefPulserMasked;
7390f655 131 if (fRefCETmean) delete fRefCETmean;
132 if (fRefCETrms) delete fRefCETrms;
133 if (fRefCEQmean) delete fRefCEQmean;
949d8707 134 if (fRefCEMasked) delete fRefCEMasked;
135 if (fRefALTROFPED) delete fRefALTROFPED;
136 if (fRefALTROZsThr) delete fRefALTROZsThr;
137 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
138 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
7390f655 139 if (fRefALTROMasked) delete fRefALTROMasked;
140 if (fRefCalibRaw) delete fRefCalibRaw;
949d8707 141 if (fCurrentRefMap) delete fCurrentRefMap;
892226be 142}
143//_____________________________________________________________________________________
144void AliTPCcalibDButil::UpdateFromCalibDB()
145{
146 //
147 // Update pointers from calibDB
148 //
1e722a63 149 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
56ce896d 150 fCalibDB->UpdateNonRec(); // load all infromation now
892226be 151 fPadNoise=fCalibDB->GetPadNoise();
152 fPedestals=fCalibDB->GetPedestals();
153 fPulserTmean=fCalibDB->GetPulserTmean();
154 fPulserTrms=fCalibDB->GetPulserTrms();
155 fPulserQmean=fCalibDB->GetPulserQmean();
156 fCETmean=fCalibDB->GetCETmean();
157 fCETrms=fCalibDB->GetCETrms();
158 fCEQmean=fCalibDB->GetCEQmean();
159 fALTROMasked=fCalibDB->GetALTROMasked();
160 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
6e7d7dc4 161 fCalibRaw=fCalibDB->GetCalibRaw();
949d8707 162 fDataQA=fCalibDB->GetDataQA();
892226be 163 UpdatePulserOutlierMap();
22c558aa 164// SetReferenceRun();
165// UpdateRefDataFromOCDB();
892226be 166}
167//_____________________________________________________________________________________
7390f655 168void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
8166d5d7 169 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
892226be 170{
171 //
172 // Process the CE data for this run
173 // the return TVectorD arrays contian the results of the fit
174 // noutliersCE contains the number of pads marked as outliers,
175 // not including masked and edge pads
176 //
177
178 //retrieve CE and ALTRO data
179 if (!fCETmean){
180 TString fitString(fitFormula);
181 fitString.ReplaceAll("++","#");
182 Int_t ndim=fitString.CountChar('#')+2;
183 fitResultsA.ResizeTo(ndim);
184 fitResultsC.ResizeTo(ndim);
185 fitResultsA.Zero();
186 fitResultsC.Zero();
187 noutliersCE=-1;
188 return;
189 }
190 noutliersCE=0;
191 //create outlier map
7390f655 192 AliTPCCalPad *out=0;
193 if (outCE) out=outCE;
194 else out=new AliTPCCalPad("outCE","outCE");
892226be 195 AliTPCCalROC *rocMasked=0x0;
196 //loop over all channels
197 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
198 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
199 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
7390f655 200 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
6e7d7dc4 201 if (!rocData) {
202 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
7390f655 203 rocOut->Add(1.);
6e7d7dc4 204 continue;
205 }
892226be 206 //add time offset to IROCs
207 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
208 rocData->Add(fIrocTimeOffset);
209 //select outliers
210 UInt_t nrows=rocData->GetNrows();
211 for (UInt_t irow=0;irow<nrows;++irow){
212 UInt_t npads=rocData->GetNPads(irow);
213 for (UInt_t ipad=0;ipad<npads;++ipad){
7390f655 214 rocOut->SetValue(irow,ipad,0);
892226be 215 //exclude masked pads
216 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
7390f655 217 rocOut->SetValue(irow,ipad,1);
892226be 218 continue;
219 }
732e90a8 220 //exclude first two rows in IROC and last two rows in OROC
221 if (iroc<36){
222 if (irow<2) rocOut->SetValue(irow,ipad,1);
223 } else {
224 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
225 }
892226be 226 //exclude edge pads
7390f655 227 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
228 Float_t valTmean=rocData->GetValue(irow,ipad);
892226be 229 //exclude values that are exactly 0
8166d5d7 230 if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
7390f655 231 rocOut->SetValue(irow,ipad,1);
892226be 232 ++noutliersCE;
233 }
234 // exclude channels with too large variations
235 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
7390f655 236 rocOut->SetValue(irow,ipad,1);
892226be 237 ++noutliersCE;
238 }
239 }
240 }
241 }
242 //perform fit
243 TMatrixD dummy;
2cb269df 244 Float_t chi2Af,chi2Cf;
245 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
246 chi2A=chi2Af;
247 chi2C=chi2Cf;
7390f655 248 if (!outCE) delete out;
892226be 249}
250//_____________________________________________________________________________________
251void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
252 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
253 Float_t &driftTimeA, Float_t &driftTimeC )
254{
255 //
256 // Calculate statistical information from the CE graphs for drift time and charge
257 //
258
259 //reset arrays
260 vecTEntries.ResizeTo(72);
261 vecTMean.ResizeTo(72);
262 vecTRMS.ResizeTo(72);
263 vecTMedian.ResizeTo(72);
264 vecQEntries.ResizeTo(72);
265 vecQMean.ResizeTo(72);
266 vecQRMS.ResizeTo(72);
267 vecQMedian.ResizeTo(72);
268 vecTEntries.Zero();
269 vecTMean.Zero();
270 vecTRMS.Zero();
271 vecTMedian.Zero();
272 vecQEntries.Zero();
273 vecQMean.Zero();
274 vecQRMS.Zero();
275 vecQMedian.Zero();
276 driftTimeA=0;
277 driftTimeC=0;
278 TObjArray *arrT=fCalibDB->GetCErocTtime();
279 TObjArray *arrQ=fCalibDB->GetCErocQtime();
280 if (arrT){
281 for (Int_t isec=0;isec<74;++isec){
282 TGraph *gr=(TGraph*)arrT->At(isec);
283 if (!gr) continue;
284 TVectorD values;
285 Int_t npoints = gr->GetN();
286 values.ResizeTo(npoints);
287 Int_t nused =0;
6e7d7dc4 288 //skip first points, theres always some problems with finding the CE position
289 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
290 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
892226be 291 values[nused]=gr->GetY()[ipoint];
292 nused++;
293 }
294 }
295 //
296 if (isec<72) vecTEntries[isec]= nused;
297 if (nused>1){
298 if (isec<72){
299 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
300 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
301 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
302 } else if (isec==72){
303 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
304 } else if (isec==73){
305 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
306 }
307 }
308 }
309 }
310 if (arrQ){
6e7d7dc4 311 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
892226be 312 TGraph *gr=(TGraph*)arrQ->At(isec);
313 if (!gr) continue;
314 TVectorD values;
315 Int_t npoints = gr->GetN();
316 values.ResizeTo(npoints);
317 Int_t nused =0;
318 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
817766d5 319 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
892226be 320 values[nused]=gr->GetY()[ipoint];
321 nused++;
322 }
323 }
324 //
6e7d7dc4 325 vecQEntries[isec]= nused;
892226be 326 if (nused>1){
327 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
328 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
329 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
330 }
331 }
332 }
333}
334
335//_____________________________________________________________________________________
336void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
337 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
93425263 338 Int_t &nonMaskedZero, Int_t &nNaN)
892226be 339{
340 //
341 // process noise data
342 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
343 // OROCs small pads [2] and OROCs large pads [3]
344 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
345 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
346 //
347
348 //set proper size and reset
349 const UInt_t infoSize=4;
350 vNoiseMean.ResizeTo(infoSize);
351 vNoiseMeanSenRegions.ResizeTo(infoSize);
352 vNoiseRMS.ResizeTo(infoSize);
353 vNoiseRMSSenRegions.ResizeTo(infoSize);
354 vNoiseMean.Zero();
355 vNoiseMeanSenRegions.Zero();
356 vNoiseRMS.Zero();
357 vNoiseRMSSenRegions.Zero();
358 nonMaskedZero=0;
93425263 359 nNaN=0;
892226be 360 //counters
361 TVectorD c(infoSize);
362 TVectorD cs(infoSize);
363 //tpc parameters
364 AliTPCParam par;
365 par.Update();
366 //retrieve noise and ALTRO data
367 if (!fPadNoise) return;
368 AliTPCCalROC *rocMasked=0x0;
369 //create IROC, OROC1, OROC2 and sensitive region masks
370 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
371 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
372 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
373 UInt_t nrows=noiseROC->GetNrows();
374 for (UInt_t irow=0;irow<nrows;++irow){
375 UInt_t npads=noiseROC->GetNPads(irow);
376 for (UInt_t ipad=0;ipad<npads;++ipad){
377 //don't use masked channels;
378 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
379 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
380 //check if noise==0
8166d5d7 381 if (noiseVal<kAlmost0) {
892226be 382 ++nonMaskedZero;
383 continue;
384 }
385 //check for nan
386 if ( !(noiseVal<10000000) ){
ef1807ef 387 AliInfo(Form("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal));
93425263 388 ++nNaN;
892226be 389 continue;
390 }
391 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
392 Int_t masksen=1; // sensitive pards are not masked (0)
393 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
394 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
395 //IROCs
396 if (irow>19&&irow<46){
397 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
398 }
399 Int_t type=1;
400 vNoiseMean[type]+=noiseVal;
401 vNoiseRMS[type]+=noiseVal*noiseVal;
402 ++c[type];
403 if (!masksen){
404 vNoiseMeanSenRegions[type]+=noiseVal;
405 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
406 ++cs[type];
407 }
408 } else {
409 //OROCs
410 //define sensive regions
411 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
412 if ( irow>75 ){
413 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
414 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
415 }
416 if ((Int_t)irow<par.GetNRowUp1()){
417 //OROC1
418 Int_t type=2;
419 vNoiseMean[type]+=noiseVal;
420 vNoiseRMS[type]+=noiseVal*noiseVal;
421 ++c[type];
422 if (!masksen){
423 vNoiseMeanSenRegions[type]+=noiseVal;
424 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
425 ++cs[type];
426 }
427 }else{
428 //OROC2
429 Int_t type=3;
430 vNoiseMean[type]+=noiseVal;
431 vNoiseRMS[type]+=noiseVal*noiseVal;
432 ++c[type];
433 if (!masksen){
434 vNoiseMeanSenRegions[type]+=noiseVal;
435 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
436 ++cs[type];
437 }
438 }
439 }
440 //whole tpc
441 Int_t type=0;
442 vNoiseMean[type]+=noiseVal;
443 vNoiseRMS[type]+=noiseVal*noiseVal;
444 ++c[type];
445 if (!masksen){
446 vNoiseMeanSenRegions[type]+=noiseVal;
447 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
448 ++cs[type];
449 }
450 }//end loop pads
451 }//end loop rows
452 }//end loop sectors (rocs)
453
454 //calculate mean and RMS
455 const Double_t verySmall=0.0000000001;
456 for (UInt_t i=0;i<infoSize;++i){
457 Double_t mean=0;
458 Double_t rms=0;
459 Double_t meanSen=0;
460 Double_t rmsSen=0;
461
462 if (c[i]>verySmall){
ef1807ef 463 AliInfo(Form("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]));
892226be 464 mean=vNoiseMean[i]/c[i];
465 rms=vNoiseRMS[i];
466 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
467 }
468 vNoiseMean[i]=mean;
469 vNoiseRMS[i]=rms;
470
471 if (cs[i]>verySmall){
472 meanSen=vNoiseMeanSenRegions[i]/cs[i];
473 rmsSen=vNoiseRMSSenRegions[i];
474 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
475 }
476 vNoiseMeanSenRegions[i]=meanSen;
477 vNoiseRMSSenRegions[i]=rmsSen;
478 }
479}
480
abb20887 481//_____________________________________________________________________________________
482void AliTPCcalibDButil::ProcessQAData(TVectorD &vQaOcc, TVectorD &vQaQtot,
483 TVectorD &vQaQmax)
484{
485 //
486 // process QA data
487 //
488 // vQaOcc/Qtot/Qmax contains the Mean occupancy/Qtot/Qmax for each sector
489 //
490
491
492 const UInt_t infoSize = 72;
493 //reset counters to error number
494 vQaOcc.ResizeTo(infoSize);
495 vQaOcc.Zero();
496 vQaQtot.ResizeTo(infoSize);
497 vQaQtot.Zero();
498 vQaQmax.ResizeTo(infoSize);
499 vQaQmax.Zero();
500 //counter
501 //retrieve pulser and ALTRO data
502
503 if (!fDataQA) {
504
505 AliInfo("No QA data");
506 return;
507 }
508 if (fDataQA->GetEventCounter()<=0) {
509
510 AliInfo("No QA data");
511 return; // no data processed
512 }
513 //
514 fDataQA->Analyse();
515
516 TVectorD normOcc(infoSize);
517 TVectorD normQ(infoSize);
518
519 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
520
ef1807ef 521 AliInfo(Form("Sector %d\n", isec));
abb20887 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
ef1807ef 533 AliInfo(Form("Nchannels %d\n", nchannels));
abb20887 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
bb22b582 837 Float_t deviation=TMath::Abs(pqRef)>1e-20?TMath::Abs(pq/pqRef-1):-999;
7390f655 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 {
668e6ebb 1231 fRefDataQA=(AliTPCdataQA*)fRefDataQA->Clone();
949d8707 1232 }
1233 delete entry;
1234 }
1235 }
1236
1237
1238//update current reference maps
1239 if (hasAnyChanged){
1240 if (fCurrentRefMap) delete fCurrentRefMap;
1241 fCurrentRefMap=(TMap*)fRefMap->Clone();
1242 }
1243}
1244//_____________________________________________________________________________________
1245AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1246{
1247 //
1248 // TObjArray object type case
1249 // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1250 //
1251 AliTPCCalPad *pad=0x0;
1252 TObjArray *arr=(TObjArray*)entry->GetObject();
1253 if (!arr){
a980538f 1254 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1255 return pad;
1256 }
1257 pad=(AliTPCCalPad*)arr->FindObject(objName);
1258 if (!pad) {
a980538f 1259 AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
949d8707 1260 return pad;
1261 }
1262 return (AliTPCCalPad*)pad->Clone();
1263}
1264//_____________________________________________________________________________________
1265AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1266{
1267 //
1268 // AliTPCCalPad object type case
1269 // cast object to a calPad and store it in 'pad'
1270 //
1271 AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1272 if (!pad) {
a980538f 1273 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1274 return 0x0;
1275 }
1276 pad=(AliTPCCalPad*)pad->Clone();
1277 return pad;
1278}
1279//_____________________________________________________________________________________
1280AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1281{
1282 //
1283 // set altro masked channel map for 'cdbPath'
1284 //
1285 AliTPCCalPad* pad=0x0;
1286 const Int_t run=GetReferenceRun(cdbPath);
1287 if (run<0) {
a980538f 1288 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1289 return pad;
1290 }
1291 AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1292 if (!entry) {
a980538f 1293 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1294 return pad;
1295 }
1296 pad=GetRefCalPad(entry,"Masked");
1297 if (pad) pad->SetNameTitle(name,name);
1298 entry->SetOwner(kTRUE);
1299 delete entry;
1300 return pad;
1301}
1302//_____________________________________________________________________________________
1303void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1304 //
1305 // Get Reference map
1306 //
1307 if (run<0) run=fCalibDB->GetRun();
1308 TString cdbPath="TPC/Calib/Ref";
1309 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1310 if (!entry) {
a980538f 1311 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
949d8707 1312 fRefMap=0;
1313 return;
1314 }
1315 entry->SetOwner(kTRUE);
1316 fRefMap=(TMap*)(entry->GetObject());
1317 AliCDBId &id=entry->GetId();
1318 fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1319}
1320//_____________________________________________________________________________________
1321Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1322{
1323 //
1324 // check whether a reference cdb entry has changed
1325 //
1326 if (!fCurrentRefMap) return kTRUE;
1327 if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1328 return kFALSE;
1329}
1330//_____________________________________________________________________________________
1331AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1332{
1333 //
1334 // get the reference AliCDBEntry for 'cdbPath'
1335 //
1336 const Int_t run=GetReferenceRun(cdbPath);
1337 if (run<0) {
a980538f 1338 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1339 return 0;
1340 }
1341 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1342 if (!entry) {
a980538f 1343 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1344 return 0;
1345 }
1346 return entry;
1347}
1348//_____________________________________________________________________________________
26d8859c 1349Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
949d8707 1350 //
1351 // Get reference run number for the specified OCDB path
1352 //
1353 if (!fCurrentRefMap) return -2;
1354 TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1355 if (!str) return -2;
58affc42 1356 return (Int_t)str->GetString().Atoi();
949d8707 1357}
1358//_____________________________________________________________________________________
26d8859c 1359Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
949d8707 1360 //
1361 // Get reference run number for the specified OCDB path
1362 //
1363 if (!fRefMap) return -1;
1364 TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1365 if (!str) return -1;
58affc42 1366 return (Int_t)str->GetString().Atoi();
949d8707 1367}
1368//_____________________________________________________________________________________
8166d5d7 1369AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
2cb269df 1370 //
1371 // Author: marian.ivanov@cern.ch
1372 //
1373 // Create outlier map for CE study
1374 // Parameters:
1375 // Return value - outlyer map
1376 // noutlyersCE - number of outlyers
1377 // minSignal - minimal total Q signal
1378 // cutRMSMin - minimal width of the signal in respect to the median
1379 // cutRMSMax - maximal width of the signal in respect to the median
1380 // cutMaxDistT - maximal deviation from time median per chamber
1381 //
1382 // Outlyers criteria:
1383 // 0. Exclude masked pads
1384 // 1. Exclude first two rows in IROC and last two rows in OROC
1385 // 2. Exclude edge pads
1386 // 3. Exclude channels with too large variations
1387 // 4. Exclude pads with too small signal
1388 // 5. Exclude signal with outlyers RMS
1389 // 6. Exclude channels to far from the chamber median
1390 noutliersCE=0;
1391 //create outlier map
1392 AliTPCCalPad *out=ceOut;
1393 if (!out) out= new AliTPCCalPad("outCE","outCE");
1394 AliTPCCalROC *rocMasked=0x0;
1395 if (!fCETmean) return 0;
1396 if (!fCETrms) return 0;
1397 if (!fCEQmean) return 0;
1398 //
1399 //loop over all channels
1400 //
1401 Double_t rmsMedian = fCETrms->GetMedian();
1402 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1403 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
9f98a33d 1404 if (!rocData) continue;
2cb269df 1405 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1406 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1407 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
1408 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1409 Double_t trocMedian = rocData->GetMedian();
1410 //
768e46f1 1411 if (!rocData || !rocCEQ || !rocCETrms || !rocData) {
2cb269df 1412 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1413 rocOut->Add(1.);
1414 continue;
1415 }
1416 //
1417 //select outliers
1418 UInt_t nrows=rocData->GetNrows();
1419 for (UInt_t irow=0;irow<nrows;++irow){
1420 UInt_t npads=rocData->GetNPads(irow);
1421 for (UInt_t ipad=0;ipad<npads;++ipad){
1422 rocOut->SetValue(irow,ipad,0);
1423 Float_t valTmean=rocData->GetValue(irow,ipad);
1424 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1425 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1426 //0. exclude masked pads
1427 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1428 rocOut->SetValue(irow,ipad,1);
1429 continue;
1430 }
1431 //1. exclude first two rows in IROC and last two rows in OROC
1432 if (iroc<36){
1433 if (irow<2) rocOut->SetValue(irow,ipad,1);
1434 } else {
1435 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1436 }
1437 //2. exclude edge pads
1438 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1439 //exclude values that are exactly 0
8166d5d7 1440 if ( TMath::Abs(valTmean)<kAlmost0) {
2cb269df 1441 rocOut->SetValue(irow,ipad,1);
1442 ++noutliersCE;
1443 }
1444 //3. exclude channels with too large variations
1445 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1446 rocOut->SetValue(irow,ipad,1);
1447 ++noutliersCE;
1448 }
1449 //
1450 //4. exclude channels with too small signal
1451 if (valQmean<minSignal) {
1452 rocOut->SetValue(irow,ipad,1);
1453 ++noutliersCE;
1454 }
1455 //
1456 //5. exclude channels with too small rms
1457 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1458 rocOut->SetValue(irow,ipad,1);
1459 ++noutliersCE;
1460 }
1461 //
1462 //6. exclude channels to far from the chamber median
1463 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1464 rocOut->SetValue(irow,ipad,1);
1465 ++noutliersCE;
1466 }
1467 }
1468 }
1469 }
1470 //
1471 return out;
1472}
1473
1474
8166d5d7 1475AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
2cb269df 1476 //
1477 // Author: marian.ivanov@cern.ch
1478 //
1479 // Create outlier map for Pulser
1480 // Parameters:
1481 // Return value - outlyer map
1482 // noutlyersPulser - number of outlyers
1483 // cutTime - absolute cut - distance to the median of chamber
1484 // cutnRMSQ - nsigma cut from median q distribution per chamber
1485 // cutnRMSrms - nsigma cut from median rms distribution
1486 // Outlyers criteria:
1487 // 0. Exclude masked pads
1488 // 1. Exclude time outlyers (default 3 time bins)
1489 // 2. Exclude q outlyers (default 5 sigma)
1490 // 3. Exclude rms outlyers (default 5 sigma)
1491 noutliersPulser=0;
1492 AliTPCCalPad *out=pulserOut;
1493 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1494 AliTPCCalROC *rocMasked=0x0;
1495 if (!fPulserTmean) return 0;
1496 if (!fPulserTrms) return 0;
1497 if (!fPulserQmean) return 0;
1498 //
1499 //loop over all channels
1500 //
1501 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1502 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1503 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1504 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1505 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1506 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1507 //
1508 Double_t rocMedianT = rocData->GetMedian();
1509 Double_t rocMedianQ = rocPulserQ->GetMedian();
1510 Double_t rocRMSQ = rocPulserQ->GetRMS();
1511 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1512 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1513 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1514 rocOut->SetValue(ichannel,0);
1515 Float_t valTmean=rocData->GetValue(ichannel);
1516 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1517 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
f03e487b 1518 Float_t valMasked =0;
1519 if (rocMasked) valMasked = rocMasked->GetValue(ichannel);
2cb269df 1520 Int_t isOut=0;
f03e487b 1521 if (valMasked>0.5) isOut=1;
2cb269df 1522 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1523 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1524 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1525 rocOut->SetValue(ichannel,isOut);
1526 if (isOut) noutliersPulser++;
1527 }
1528 }
1529 return out;
1530}
1531
1532
1533AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1534 //
1535 // Author : Marian Ivanov
1536 // Create pad time0 correction map using information from the CE and from pulser
1537 //
1538 //
1539 // Return PadTime0 to be used for time0 relative alignment
1540 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1541 // using $ALICE_ROOT/TPC/script/gui application
1542 //
1543 // fitResultsA - fitParameters A side
1544 // fitResultsC - fitParameters C side
1545 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1546 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1547 //
1548 //
1549 // Algorithm:
1550 // 1. Find outlier map for CE
1551 // 2. Find outlier map for Pulser
1552 // 3. Replace outlier by median at given sector (median without outliers)
1553 // 4. Substract from the CE data pulser
1554 // 5. Fit the CE with formula
1555 // 5.1) (IROC-OROC) offset
1556 // 5.2) gx
1557 // 5.3) gy
1558 // 5.4) (lx-xmid)
1559 // 5.5) (IROC-OROC)*(lx-xmid)
1560 // 5.6) (ly/lx)^2
1561 // 6. Substract gy fit dependence from the CE data
1562 // 7. Add pulser back to CE data
1563 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1564 // 9. return CE data
1565 //
1566 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1567 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1568
1569 // fit formula
1570 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)";
1571 // output for fit formula
1572 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)";
1573 // gy part of formula
1574 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)";
1575 //
1576 //
1577 if (!fCETmean) return 0;
1578 Double_t pgya,pgyc,pchi2a,pchi2c;
1579 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1580 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1581
1582 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1583 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1584 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1585 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1586 padPulser->SetName("padPulser");
1587 padPulserOut->SetName("padPulserOut");
1588 padCE->SetName("padCE");
1589 padCEIn->SetName("padCEIn");
1590 padCEOut->SetName("padCEOut");
1591 padOut->SetName("padOut");
1592
1593 //
1594 // make combined outlyers map
1595 // and replace outlyers in maps with median for chamber
1596 //
1597 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1598 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1599 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1600 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1601 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1602 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1603 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1604 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1605 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1606 if (rocPulserOut->GetValue(ichannel)>0) {
1607 rocPulser->SetValue(ichannel,pulserMedian);
1608 rocOut->SetValue(ichannel,1);
1609 }
1610 if (rocCEOut->GetValue(ichannel)>0) {
1611 rocCE->SetValue(ichannel,ceMedian);
1612 rocOut->SetValue(ichannel,1);
1613 }
1614 }
1615 }
1616 //
1617 // remove pulser time 0
1618 //
1619 padCE->Add(padPulser,-1);
1620 //
1621 // Make fits
1622 //
1623 TMatrixD dummy;
1624 Float_t chi2Af,chi2Cf;
1625 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1626 chi2A=chi2Af;
1627 chi2C=chi2Cf;
1628 //
1629 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1630 padCEFitGY->SetName("padCEFitGy");
1631 //
1632 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1633 padCEFit->SetName("padCEFit");
1634 //
1635 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1636 padCEDiff->SetName("padCEDiff");
1637 padCEDiff->Add(padCEFit,-1.);
1638 //
1639 //
1640 padCE->Add(padCEFitGY,-1.);
1641
1642 padCE->Add(padPulser,1.);
1643 Double_t padmedian = padCE->GetMedian();
1644 padCE->Add(-padmedian); // normalize to median
1645 //
1646 // Replace outliers by fit value - median of diff per given chamber -GY fit
1647 //
1648 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1649 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1650 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1651 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1652 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1653 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1654 //
1655 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1656 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1657 if (rocOut->GetValue(ichannel)==0) continue;
1658 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1659 rocCE->SetValue(ichannel,value);
1660 }
1661 }
1662 //
1663 //
1664 if (dumpfile){
1665 //dump to the file - result can be visualized
1666 AliTPCPreprocessorOnline preprocesor;
1667 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1668 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1669 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1670 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1671 //
1672 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1673 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1674 //
1675 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1676 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1677 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1678 preprocesor.DumpToFile(dumpfile);
1679 }
1680 delete padPulser;
1681 delete padPulserOut;
1682 delete padCEIn;
1683 delete padCEOut;
1684 delete padOut;
1685 delete padCEDiff;
1686 delete padCEFitGY;
1687 return padCE;
1688}
1689
817766d5 1690
1691
1692
1693
1694Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1695 //
1696 // find the closest point to xref in x direction
1697 // return dx and value
5647625c 1698 dx = 0;
1699 y = 0;
1700
1701 if(!graph) return 0;
1702 if(graph->GetN() < 1) return 0;
1703
817766d5 1704 Int_t index=0;
1705 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1706 if (index<0) index=0;
5647625c 1707 if(graph->GetN()==1) {
1708 dx = xref-graph->GetX()[index];
1709 }
1710 else {
1711 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1712 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1713 dx = xref-graph->GetX()[index];
1714 }
817766d5 1715 y = graph->GetY()[index];
1716 return index;
1717}
1718
817766d5 1719Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1720 //
1721 // Get the correction of the trigger offset
1722 // combining information from the laser track calibration
1723 // and from cosmic calibration
1724 //
1725 // run - run number
1726 // timeStamp - tim stamp in seconds
1727 // deltaT - integration period to calculate offset
1728 // deltaTLaser -max validity of laser data
1729 // valType - 0 - median, 1- mean
1730 //
1731 // Integration vaues are just recomendation - if not possible to get points
1732 // automatically increase the validity by factor 2
1733 // (recursive algorithm until one month of data taking)
1734 //
1735 //
1736 const Float_t kLaserCut=0.0005;
aec07590 1737 const Int_t kMaxPeriod=3600*24*30*12; // one year max
817766d5 1738 const Int_t kMinPoints=20;
1739 //
1740 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1741 if (!array) {
1742 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1743 }
1744 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1745 if (!array) return 0;
1746 //
1747 TGraphErrors *laserA[3]={0,0,0};
1748 TGraphErrors *laserC[3]={0,0,0};
1749 TGraphErrors *cosmicAll=0;
1750 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1751 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1752 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1753 //
1754 //
1755 if (!cosmicAll) return 0;
1756 Int_t nmeasC=cosmicAll->GetN();
1757 Float_t *tdelta = new Float_t[nmeasC];
1758 Int_t nused=0;
1759 for (Int_t i=0;i<nmeasC;i++){
1760 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1761 Float_t ccosmic=cosmicAll->GetY()[i];
1762 Double_t yA=0,yC=0,dA=0,dC=0;
1763 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1764 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1765 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1766 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1767 //
1768 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1769 Float_t claser=0;
1770 if (TMath::Abs(yA-yC)<kLaserCut) {
1771 claser=(yA-yC)*0.5;
1772 }else{
1773 if (i%2==0) claser=yA;
1774 if (i%2==1) claser=yC;
1775 }
1776 tdelta[nused]=ccosmic-claser;
1777 nused++;
1778 }
768e46f1 1779 if (nused<kMinPoints &&deltaT<kMaxPeriod) {
1780 delete [] tdelta;
1781 return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1782 }
aec07590 1783 if (nused<kMinPoints) {
f03e487b 1784 delete [] tdelta;
6157d837 1785 //AliWarning("AliFatal: No time offset calibration available\n");
aec07590 1786 return 0;
1787 }
817766d5 1788 Double_t median = TMath::Median(nused,tdelta);
1789 Double_t mean = TMath::Mean(nused,tdelta);
a9f487a1 1790 delete [] tdelta;
817766d5 1791 return (valType==0) ? median:mean;
1792}
1793
a23ba1c3 1794Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
817766d5 1795 //
1796 // Get the correction of the drift velocity
1797 // combining information from the laser track calibration
1798 // and from cosmic calibration
1799 //
a23ba1c3 1800 // dist - return value - distance to closest point in graph
817766d5 1801 // run - run number
1802 // timeStamp - tim stamp in seconds
1803 // deltaT - integration period to calculate time0 offset
1804 // deltaTLaser -max validity of laser data
1805 // valType - 0 - median, 1- mean
1806 //
1807 // Integration vaues are just recomendation - if not possible to get points
1808 // automatically increase the validity by factor 2
1809 // (recursive algorithm until one month of data taking)
1810 //
1811 //
1812 //
1813 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1814 if (!array) {
1815 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1816 }
1817 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1818 if (!array) return 0;
1819 TGraphErrors *cosmicAll=0;
1820 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1821 if (!cosmicAll) return 0;
a23ba1c3 1822 Double_t grY=0;
1823 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1824
817766d5 1825 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
5647625c 1826 Double_t vcosmic = AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1e722a63 1827 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
817766d5 1828 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
cc65e4f5 1829 return vcosmic-t0;
817766d5 1830
1831 /*
1832 Example usage:
1833
1834 Int_t run=89000
1835 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1836 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1837 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1838 //
1839 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1840 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1841 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1842 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1843
1844 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1845 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1846
1847 */
1848
1849}
1850
949d8707 1851const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1852{
1853 //
1854 // Create a default name for the gui file
1855 //
1856
1857 return Form("guiRefTreeRun%s.root",GetRefValidity());
1858}
1859
1860Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1861{
1862 //
1863 // Create a gui reference tree
1864 // if dirname and filename are empty default values will be used
1865 // this is the recommended way of using this function
1866 // it allows to check whether a file with the given run validity alredy exists
1867 //
1868 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1869 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1870 return kFALSE;
1871 }
1872
1873 TString file=filename;
1874 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1875
1876 AliTPCPreprocessorOnline prep;
1877 //noise and pedestals
1878 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1879 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1880 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1881 //pulser data
1882 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1883 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1884 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1885 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1886 //CE data
1887 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1888 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1889 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1890 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1891 //Altro data
1892 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1893 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1894 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1895 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1896 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1897 //QA
1898 AliTPCdataQA *dataQA=fRefDataQA;
1899 if (dataQA) {
1900 if (dataQA->GetNLocalMaxima())
1901 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1902 if (dataQA->GetMaxCharge())
1903 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1904 if (dataQA->GetMeanCharge())
1905 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1906 if (dataQA->GetNoThreshold())
1907 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1908 if (dataQA->GetNTimeBins())
1909 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1910 if (dataQA->GetNPads())
1911 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1912 if (dataQA->GetTimePosition())
1913 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1914 }
1915 prep.DumpToFile(file.Data());
1916 return kTRUE;
1917}
1918
a23ba1c3 1919Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1920 //
78f17711 1921 // Get the correction of the drift velocity using the offline laser tracks calbration
a23ba1c3 1922 //
1923 // run - run number
1924 // timeStamp - tim stamp in seconds
1925 // deltaT - integration period to calculate time0 offset
1926 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1927 // Note in case no data form both A and C side - the value from active side used
1928 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
78f17711 1929
1930 return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1931}
1932
1933Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1934 //
1935 // Get the correction of the drift velocity using the online laser tracks calbration
1936 //
1937 // run - run number
1938 // timeStamp - tim stamp in seconds
1939 // deltaT - integration period to calculate time0 offset
1940 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1941 // Note in case no data form both A and C side - the value from active side used
1942 TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1943
1944 Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1945 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1946 if (!param) return 0;
1947
1948 //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1949 dv*=param->GetDriftV()/2.61301900000000000e+06;
1950 if (dv>1e-20) dv=1/dv-1;
1951 else return 0;
1952 // T/P correction
1953 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1954
1955 AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
2cd3a697 1956 AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
78f17711 1957
1958 Double_t corrPTA=0;
1959 Double_t corrPTC=0;
1960
1961 if (temp&&press) {
1962 AliTPCCalibVdrift corr(temp,press,0);
1963 corrPTA=corr.GetPTRelative(timeStamp,0);
1964 corrPTC=corr.GetPTRelative(timeStamp,1);
1965 }
1966
1967 if (side==0) dv -= corrPTA;
1968 if (side==1) dv -= corrPTC;
1969 if (side==2) dv -= (corrPTA+corrPTC)/2;
1970
1971 return dv;
1972}
1973
1974Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1975 Int_t side, TObjArray * const array){
1976 //
1977 // common drift velocity retrieval for online and offline method
1978 //
a23ba1c3 1979 TGraphErrors *grlaserA=0;
1980 TGraphErrors *grlaserC=0;
1981 Double_t vlaserA=0, vlaserC=0;
1982 if (!array) return 0;
1983 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1984 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1985 Double_t deltaY;
5647625c 1986 if (grlaserA && grlaserA->GetN()>0) {
a23ba1c3 1987 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1988 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1989 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1990 }
5647625c 1991 if (grlaserC && grlaserC->GetN()>0) {
a23ba1c3 1992 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1993 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1994 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1995 }
1996 if (side==0) return vlaserA;
1997 if (side==1) return vlaserC;
1998 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1999 if (!grlaserA) return vlaserC;
2000 if (!grlaserC) return vlaserA;
2001 return mdrift;
2002}
2003
2004
2005Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
2006 //
2007 // Get the correction of the drift velocity using the CE laser data
2008 // combining information from the CE, laser track calibration
2009 // and P/T calibration
2010 //
2011 // run - run number
2012 // timeStamp - tim stamp in seconds
2013 // deltaT - integration period to calculate time0 offset
2014 // side - 0 - A side, 1 - C side, 2 - mean from both sides
2015 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
0d1b4cf8 2016 if (!arrT) return 0;
a23ba1c3 2017 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
2018 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2019 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2020 //
2021 //
2022 Double_t corrPTA = 0, corrPTC=0;
2023 Double_t ltime0A = 0, ltime0C=0;
2024 Double_t gry=0;
2025 Double_t corrA=0, corrC=0;
2026 Double_t timeA=0, timeC=0;
5647625c 2027 const Double_t kEpsilon = 0.00001;
a23ba1c3 2028 TGraph *graphA = (TGraph*)arrT->At(72);
2029 TGraph *graphC = (TGraph*)arrT->At(73);
2030 if (!graphA && !graphC) return 0.;
2031 if (graphA &&graphA->GetN()>0) {
2032 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2033 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2034 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2035 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2036 if(ltime0A < kEpsilon) return 0;
a23ba1c3 2037 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2038 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2039 corrA-=corrPTA;
2040 }
2041 if (graphC&&graphC->GetN()>0){
2042 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2043 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2044 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2045 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2046 if(ltime0C < kEpsilon) return 0;
2047if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2048 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2049 corrC-=corrPTC;
2050 }
2051
2052 if (side ==0 ) return corrA;
2053 if (side ==1 ) return corrC;
2054 Double_t corrM= (corrA+corrC)*0.5;
2055 if (!graphA) corrM=corrC;
2056 if (!graphC) corrM=corrA;
2057 return corrM;
2058}
2059
cc65e4f5 2060Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2061 //
2062 // return drift velocity using the TPC-ITS matchin method
2063 // return also distance to the closest point
2064 //
2065 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2066 TGraphErrors *graph=0;
2067 dist=0;
2068 if (!array) return 0;
5647625c 2069 //array->ls();
cc65e4f5 2070 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
3760239f 2071 if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_TOFB_TPC_DRIFTVD");
cc65e4f5 2072 if (!graph) return 0;
2073 Double_t deltaY;
2074 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2075 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2076 return value;
2077}
2078
2079Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2080 //
2081 // Get time dependent time 0 (trigger delay in cm) correction
2082 // Arguments:
2083 // timestamp - timestamp
2084 // run - run number
2085 //
2086 // Notice - Extrapolation outside of calibration range - using constant function
2087 //
2088 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2089 TGraphErrors *graph=0;
2090 dist=0;
2091 if (!array) return 0;
3760239f 2092 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
2093 if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
cc65e4f5 2094 if (!graph) return 0;
2095 Double_t deltaY;
2096 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2097 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2098 return value;
2099}
2100
2101
a23ba1c3 2102
2103
2104
2105Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2106 //
2107 // VERY obscure method - we need something in framework
2108 // Find the TPC runs with temperature OCDB entry
2109 // cache the start and end of the run
2110 //
2111 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2112 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2113 if (!storage) return 0;
2114 TString path=storage->GetURI();
2115 TString runsT;
2116 {
2117 TString command;
2118 if (path.Contains("local")){ // find the list if local system
2119 path.ReplaceAll("local://","");
2120 path+="TPC/Calib/Temperature";
2121 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2122 }
2123 runsT=gSystem->GetFromPipe(command);
2124 }
2125 TObjArray *arr= runsT.Tokenize("r");
2126 if (!arr) return 0;
2127 //
2128 TArrayI indexes(arr->GetEntries());
2129 TArrayI runs(arr->GetEntries());
2130 Int_t naccept=0;
2131 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2132 Int_t irunN = atoi(arr->At(irun)->GetName());
2133 if (irunN<startRun) continue;
2134 if (irunN>stopRun) continue;
2135 runs[naccept]=irunN;
2136 naccept++;
2137 }}
09d5920f 2138 delete arr;
a23ba1c3 2139 fRuns.Set(naccept);
2140 fRunsStart.Set(fRuns.fN);
2141 fRunsStop.Set(fRuns.fN);
2142 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2143 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2144
2145 //
2146 AliCDBEntry * entry = 0;
2147 {for (Int_t irun=0;irun<fRuns.fN; irun++){
2148 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2149 if (!entry) continue;
2150 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2151 if (!tmpRun) continue;
2152 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2153 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
ef1807ef 2154 //AliInfo(Form("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec()));
a23ba1c3 2155 }}
2156 return fRuns.fN;
2157}
2158
2159
2160Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2161 //
2162 // binary search - find the run for given time stamp
2163 //
2164 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2165 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2166 Int_t cindex = -1;
2167 for (Int_t index=index0; index<=index1; index++){
2168 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2169 if (debug) {
ef1807ef 2170 AliInfo(Form("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime));
a23ba1c3 2171 }
2172 }
2173 if (cindex<0) cindex =(index0+index1)/2;
2174 if (cindex<0) {
2175 return 0;
2176 }
2177 return fRuns[cindex];
2178}
2179
2180
2181
2182
2183
2184TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2185 //
2186 // filter outlyer measurement
2187 // Only points around median +- sigmaCut filtered
2188 if (!graph) return 0;
2189 Int_t kMinPoints=2;
2190 Int_t npoints0 = graph->GetN();
2191 Int_t npoints=0;
2192 Float_t rmsY=0;
a23ba1c3 2193 //
2194 //
2195 if (npoints0<kMinPoints) return 0;
a9f487a1 2196
2197 Double_t *outx=new Double_t[npoints0];
2198 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2199 for (Int_t iter=0; iter<3; iter++){
2200 npoints=0;
2201 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2202 if (graph->GetY()[ipoint]==0) continue;
2203 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2204 outx[npoints] = graph->GetX()[ipoint];
2205 outy[npoints] = graph->GetY()[ipoint];
2206 npoints++;
2207 }
2208 if (npoints<=1) break;
2209 medianY =TMath::Median(npoints,outy);
2210 rmsY =TMath::RMS(npoints,outy);
2211 }
2212 TGraph *graphOut=0;
2213 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2214 delete [] outx;
2215 delete [] outy;
a23ba1c3 2216 return graphOut;
2217}
2218
2219
2220TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2221 //
2222 // filter outlyer measurement
2223 // Only points around median +- cut filtered
2224 if (!graph) return 0;
2225 Int_t kMinPoints=2;
2226 Int_t npoints0 = graph->GetN();
2227 Int_t npoints=0;
829455ad 2228// Float_t rmsY=0;
a23ba1c3 2229 //
2230 //
2231 if (npoints0<kMinPoints) return 0;
a9f487a1 2232
2233 Double_t *outx=new Double_t[npoints0];
2234 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2235 for (Int_t iter=0; iter<3; iter++){
2236 npoints=0;
2237 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2238 if (graph->GetY()[ipoint]==0) continue;
2239 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2240 outx[npoints] = graph->GetX()[ipoint];
2241 outy[npoints] = graph->GetY()[ipoint];
2242 npoints++;
2243 }
2244 if (npoints<=1) break;
2245 medianY =TMath::Median(npoints,outy);
829455ad 2246// rmsY =TMath::RMS(npoints,outy);
a23ba1c3 2247 }
2248 TGraph *graphOut=0;
2249 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2250 delete [] outx;
2251 delete [] outy;
a23ba1c3 2252 return graphOut;
2253}
2254
2255
2256
8166d5d7 2257TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
a23ba1c3 2258 //
2259 // filter outlyer measurement
2260 // Only points with normalized errors median +- sigmaCut filtered
2261 //
2262 Int_t kMinPoints=10;
2263 Int_t npoints0 = graph->GetN();
2264 Int_t npoints=0;
2265 Float_t medianErr=0, rmsErr=0;
a9f487a1 2266 //
2267 //
2268 if (npoints0<kMinPoints) return 0;
2269
a23ba1c3 2270 Double_t *outx=new Double_t[npoints0];
2271 Double_t *outy=new Double_t[npoints0];
2272 Double_t *erry=new Double_t[npoints0];
2273 Double_t *nerry=new Double_t[npoints0];
2274 Double_t *errx=new Double_t[npoints0];
a9f487a1 2275
a23ba1c3 2276 for (Int_t iter=0; iter<3; iter++){
2277 npoints=0;
2278 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2279 nerry[npoints] = graph->GetErrorY(ipoint);
2280 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2281 erry[npoints] = graph->GetErrorY(ipoint);
2282 outx[npoints] = graph->GetX()[ipoint];
2283 outy[npoints] = graph->GetY()[ipoint];
2284 errx[npoints] = graph->GetErrorY(ipoint);
2285 npoints++;
2286 }
2287 if (npoints==0) break;
2288 medianErr=TMath::Median(npoints,erry);
2289 medianY =TMath::Median(npoints,outy);
2290 rmsErr =TMath::RMS(npoints,erry);
2291 }
2292 TGraphErrors *graphOut=0;
2293 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2294 delete []outx;
2295 delete []outy;
a23ba1c3 2296 delete []erry;
a9f487a1 2297 delete []nerry;
2298 delete []errx;
a23ba1c3 2299 return graphOut;
2300}
2301
2302
2303void AliTPCcalibDButil::Sort(TGraph *graph){
2304 //
2305 // sort array - neccessay for approx
2306 //
2307 Int_t npoints = graph->GetN();
2308 Int_t *indexes=new Int_t[npoints];
2309 Double_t *outx=new Double_t[npoints];
2310 Double_t *outy=new Double_t[npoints];
2311 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2312 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2313 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2314 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2315 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2316
a9f487a1 2317 delete [] indexes;
2318 delete [] outx;
2319 delete [] outy;
a23ba1c3 2320}
2321void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2322 //
2323 // smmoth graph - mean on the interval
2324 //
2325 Sort(graph);
2326 Int_t npoints = graph->GetN();
2327 Double_t *outy=new Double_t[npoints];
0d1b4cf8 2328
a23ba1c3 2329 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2330 Double_t lx=graph->GetX()[ipoint];
2331 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2332 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2333 if (index0<0) index0=0;
2334 if (index1>=npoints-1) index1=npoints-1;
2335 if ((index1-index0)>1){
2336 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2337 }else{
2338 outy[ipoint]=graph->GetY()[ipoint];
2339 }
2340 }
0d1b4cf8 2341 // TLinearFitter fitter(3,"pol2");
2342// for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2343// Double_t lx=graph->GetX()[ipoint];
2344// Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2345// Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2346// if (index0<0) index0=0;
2347// if (index1>=npoints-1) index1=npoints-1;
2348// fitter.ClearPoints();
2349// for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2350// if ((index1-index0)>1){
2351// outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2352// }else{
2353// outy[ipoint]=graph->GetY()[ipoint];
2354// }
2355// }
2356
2357
2358
a23ba1c3 2359 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2360 graph->GetY()[ipoint] = outy[ipoint];
2361 }
2362 delete[] outy;
2363}
2364
8166d5d7 2365Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
a23ba1c3 2366 //
2367 // Use constant interpolation outside of range
2368 //
2369 if (!graph) {
ef1807ef 2370 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
a23ba1c3 2371 return 0;
2372 }
5647625c 2373
a23ba1c3 2374 if (graph->GetN()<1){
ef1807ef 2375 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
a23ba1c3 2376 return 0;
2377 }
5647625c 2378
2379
a23ba1c3 2380 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2381 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
5647625c 2382
ef1807ef 2383 // AliInfo(Form("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref)));
5647625c 2384
2385 if(graph->GetN()==1)
2386 return graph->Eval(graph->GetX()[0]);
2387
2388
2389 return graph->Eval(xref);
a23ba1c3 2390}
2391
108a67db 2392Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2393 //
2394 // Use constant interpolation outside of range also for spline fits
2395 //
2396 if (!graph) {
ef1807ef 2397 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
108a67db 2398 return 0;
2399 }
2400 if (graph->GetKnots()<1){
ef1807ef 2401 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph");
108a67db 2402 return 0;
2403 }
2404 if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2405 if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2406 return graph->Eval( xref);
2407}
2408
1e722a63 2409Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2410 //
2411 // Filter DCS sensor information
2412 // ymin - minimal value
2413 // ymax - max value
2414 // maxdy - maximal deirivative
2415 // sigmaCut - cut on values and derivative in terms of RMS distribution
2416 // Return value - accepted fraction
2417 //
2418 // Algorithm:
2419 //
2420 // 0. Calculate median and rms of values in specified range
2421 // 1. Filter out outliers - median+-sigmaCut*rms
2422 // values replaced by median
2423 //
2424 AliSplineFit * fit = sensor->GetFit();
2425 if (!fit) return 0.;
2426 Int_t nknots = fit->GetKnots();
2427 if (nknots==0) {
2428 delete fit;
2429 sensor->SetFit(0);
2430 return 0;
2431 }
2432 //
2433 Double_t *yin0 = new Double_t[nknots];
2434 Double_t *yin1 = new Double_t[nknots];
2435 Int_t naccept=0;
2436
2437 for (Int_t iknot=0; iknot< nknots; iknot++){
2438 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2439 yin0[naccept] = fit->GetY0()[iknot];
2440 yin1[naccept] = fit->GetY1()[iknot];
2441 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2442 naccept++;
2443 }
2444 }
2445 if (naccept<1) {
2446 delete fit;
2447 sensor->SetFit(0);
a9f487a1 2448 delete [] yin0;
2449 delete [] yin1;
1e722a63 2450 return 0.;
2451 }
0d1b4cf8 2452
1e722a63 2453 Double_t medianY0=0, medianY1=0;
2454 Double_t rmsY0 =0, rmsY1=0;
2455 medianY0 = TMath::Median(naccept, yin0);
2456 medianY1 = TMath::Median(naccept, yin1);
2457 rmsY0 = TMath::RMS(naccept, yin0);
2458 rmsY1 = TMath::RMS(naccept, yin1);
2459 naccept=0;
2460 //
2461 // 1. Filter out outliers - median+-sigmaCut*rms
2462 // values replaced by median
2463 // if replaced the derivative set to 0
2464 //
2465 for (Int_t iknot=0; iknot< nknots; iknot++){
2466 Bool_t isOK=kTRUE;
2467 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2468 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2469 if (nknots<2) fit->GetY1()[iknot]=0;
2470 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2471 if (!isOK){
2472 fit->GetY0()[iknot]=medianY0;
2473 fit->GetY1()[iknot]=0;
2474 }else{
2475 naccept++;
2476 }
2477 }
2478 delete [] yin0;
2479 delete [] yin1;
2480 return Float_t(naccept)/Float_t(nknots);
2481}
2482
2483Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2484 //
2485 // Filter temperature array
2486 // tempArray - array of temperatures -
2487 // ymin - minimal accepted temperature - default 15
2488 // ymax - maximal accepted temperature - default 22
2489 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2490 // return value - fraction of filtered sensors
2491 const Double_t kMaxDy=0.1;
2492 Int_t nsensors=tempArray->NumSensors();
2493 if (nsensors==0) return 0.;
2494 Int_t naccept=0;
2495 for (Int_t isensor=0; isensor<nsensors; isensor++){
2496 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2497 if (!sensor) continue;
1e722a63 2498 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2499 if (sensor->GetFit()==0){
0d1b4cf8 2500 //delete sensor;
1e722a63 2501 tempArray->RemoveSensorNum(isensor);
2502 }else{
2503 naccept++;
2504 }
2505 }
2506 return Float_t(naccept)/Float_t(nsensors);
2507}
2508
a23ba1c3 2509
8166d5d7 2510void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2511 //
2512 // Filter CE data
1e722a63 2513 // Input parameters:
2514 // deltaT - smoothing window (in seconds)
2515 // cutAbs - max distance of the time info to the median (in time bins)
2516 // cutSigma - max distance (in the RMS)
2517 // pcstream - optional debug streamer to store original and filtered info
2518 // Hardwired parameters:
2519 // kMinPoints =10; // minimal number of points to define the CE
2520 // kMinSectors=12; // minimal number of sectors to define sideCE
2521 // Algorithm:
2522 // 0. Filter almost emty graphs (kMinPoints=10)
2523 // 1. calculate median and RMS per side
2524 // 2. Filter graphs - in respect with side medians
2525 // - cutAbs and cutDelta used
2526 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2527 // 4. Calculate mean for A side and C side
2528 //
2529 const Int_t kMinPoints =10; // minimal number of points to define the CE
2530 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2531 const Int_t kMinTime =400; // minimal arrival time of CE
2532 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
a23ba1c3 2533 Double_t medianY=0;
1e722a63 2534 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
a23ba1c3 2535 if (!cearray) return;
1e722a63 2536 Double_t tmin=-1;
2537 Double_t tmax=-1;
2538 //
2539 //
a23ba1c3 2540 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2cd3a697 2541 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
a23ba1c3 2542 if ( tempMapCE && cavernPressureCE){
1e722a63 2543 //
c5dfa2dc 2544 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2545 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2546 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2547 Bool_t isOK=kTRUE;
1e722a63 2548 if (isOK) {
2549 // recalculate P/T correction map for time of the CE
2550 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2551 driftCalib->SetName("driftPTCE");
2552 driftCalib->SetTitle("driftPTCE");
2553 cearray->AddLast(driftCalib);
2554 }
a23ba1c3 2555 }
1e722a63 2556 //
2557 // 0. Filter almost emty graphs
2558 //
a23ba1c3 2559
1e722a63 2560 for (Int_t i=0; i<72;i++){
a23ba1c3 2561 TGraph *graph= (TGraph*)arrT->At(i);
29797f44 2562 if (!graph) continue;
2563 graph->Sort();
a23ba1c3 2564 if (graph->GetN()<kMinPoints){
2565 arrT->AddAt(0,i);
2566 delete graph; // delete empty graph
2567 continue;
2568 }
1e722a63 2569 if (tmin<0) tmin = graph->GetX()[0];
2570 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2571 //
2572 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2573 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2574 }
2575 //
2576 // 1. calculate median and RMS per side
2577 //
2578 TArrayF arrA(100000), arrC(100000);
2579 Int_t nA=0, nC=0;
2580 Double_t medianA=0, medianC=0;
2581 Double_t rmsA=0, rmsC=0;
2582 for (Int_t isec=0; isec<72;isec++){
2583 TGraph *graph= (TGraph*)arrT->At(isec);
2584 if (!graph) continue;
2585 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2586 if (graph->GetY()[ipoint]<kMinTime) continue;
2587 if (nA>=arrA.fN) arrA.Set(nA*2);
2588 if (nC>=arrC.fN) arrC.Set(nC*2);
2589 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2590 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2591 }
2592 }
2593 if (nA>0){
2594 medianA=TMath::Median(nA,arrA.fArray);
2595 rmsA =TMath::RMS(nA,arrA.fArray);
2596 }
2597 if (nC>0){
2598 medianC=TMath::Median(nC,arrC.fArray);
2599 rmsC =TMath::RMS(nC,arrC.fArray);
2600 }
2601 //
2602 // 2. Filter graphs - in respect with side medians
2603 //
2604 TArrayD vecX(100000), vecY(100000);
2605 for (Int_t isec=0; isec<72;isec++){
2606 TGraph *graph= (TGraph*)arrT->At(isec);
2607 if (!graph) continue;
2608 Double_t median = (isec%36<18) ? medianA: medianC;
2609 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2610 Int_t naccept=0;
29797f44 2611 // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2612 for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
1e722a63 2613 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2614 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2615 vecX[naccept]= graph->GetX()[ipoint];
2616 vecY[naccept]= graph->GetY()[ipoint];
2617 naccept++;
2618 }
2619 if (naccept<kMinPoints){
2620 arrT->AddAt(0,isec);
2621 delete graph; // delete empty graph
2622 continue;
2623 }
2624 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2625 delete graph;
2626 arrT->AddAt(graph2,isec);
2627 }
2628 //
2629 // 3. Cut in respect wit the graph median
2630 //
2631 for (Int_t i=0; i<72;i++){
2632 TGraph *graph= (TGraph*)arrT->At(i);
2633 if (!graph) continue;
2634 //
2635 // filter in range
2636 //
2637 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
a23ba1c3 2638 if (!graphTS0) continue;
2639 if (graphTS0->GetN()<kMinPoints) {
2640 delete graphTS0;
2641 delete graph;
2642 arrT->AddAt(0,i);
2643 continue;
2644 }
1e722a63 2645 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
f03e487b 2646 if (!graphTS) continue;
a23ba1c3 2647 graphTS->Sort();
2648 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2649 if (pcstream){
2650 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2651 (*pcstream)<<"filterCE"<<
2652 "run="<<run<<
2653 "isec="<<i<<
2654 "mY="<<medianY<<
2655 "graph.="<<graph<<
2656 "graphTS0.="<<graphTS0<<
2657 "graphTS.="<<graphTS<<
2658 "\n";
2659 }
2660 delete graphTS0;
a23ba1c3 2661 arrT->AddAt(graphTS,i);
2662 delete graph;
2663 }
1e722a63 2664 //
2665 // Recalculate the mean time A side C side
2666 //
2667 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2668 Int_t meanPoints=(nA+nC)/72; // mean number of points
2669 for (Int_t itime=0; itime<200; itime++){
2670 nA=0, nC=0;
2671 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2672 for (Int_t i=0; i<72;i++){
2673 TGraph *graph= (TGraph*)arrT->At(i);
2674 if (!graph) continue;
2675 if (graph->GetN()<(meanPoints/4)) continue;
2676 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2677 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2678 }
2679 xA[itime]=time;
2680 xC[itime]=time;
2681 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2682 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2683 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2684 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2685 }
2686 //
2687 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2688 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2689 if (pcstream){
2690 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2691 (*pcstream)<<"filterAC"<<
2692 "run="<<run<<
2693 "nA="<<nA<<
2694 "nC="<<nC<<
2695 "rmsTA="<<rmsTA<<
2696 "rmsTC="<<rmsTC<<
2697 "\n";
2698 }
2699 //
2700 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2701 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2702 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2703 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2704 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2705 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2706 delete grA;
2707 delete grC;
2708 if (nA<kMinSectors) arrT->AddAt(0,72);
2709 else arrT->AddAt(graphTSA,72);
2710 if (nC<kMinSectors) arrT->AddAt(0,73);
2711 else arrT->AddAt(graphTSC,73);
a23ba1c3 2712}
2713
2714
8166d5d7 2715void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2716 //
2717 // Filter Drift velocity measurement using the tracks
2718 // 0. remove outlyers - error based
2719 // cutSigma
2720 //
2721 //
2722 const Int_t kMinPoints=1; // minimal number of points to define value
2723 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2724 Double_t medianY=0;
2725 if (!arrT) return;
2726 for (Int_t i=0; i<arrT->GetEntries();i++){
14301155 2727 TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
a23ba1c3 2728 if (!graph) continue;
2729 if (graph->GetN()<kMinPoints){
2730 delete graph;
2731 arrT->AddAt(0,i);
2732 continue;
2733 }
5647625c 2734 TGraphErrors *graph2 = NULL;
2735 if(graph->GetN()<10) {
2736 graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2737 if (!graph2) {
2738 delete graph; arrT->AddAt(0,i); continue;
2739 }
2740 }
2741 else {
2742 graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2743 if (!graph2) {
2744 delete graph; arrT->AddAt(0,i); continue;
2745 }
a23ba1c3 2746 }
2747 if (graph2->GetN()<1) {
2748 delete graph; arrT->AddAt(0,i); continue;
2749 }
2750 graph2->SetName(graph->GetName());
2751 graph2->SetTitle(graph->GetTitle());
2752 arrT->AddAt(graph2,i);
2753 if (pcstream){
2754 (*pcstream)<<"filterTracks"<<
2755 "run="<<run<<
2756 "isec="<<i<<
2757 "mY="<<medianY<<
2758 "graph.="<<graph<<
2759 "graph2.="<<graph2<<
2760 "\n";
2761 }
2762 delete graph;
2763 }
2764}
2765
2766
a23ba1c3 2767
2768
2769
2770Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2771 //
2772 //
2773 // get laser time offset
2774 // median around timeStamp+-deltaT
2775 // QA - chi2 needed for later usage - to be added
2776 // - currently cut on error
2777 //
2778 Int_t kMinPoints=1;
2779 Double_t kMinDelay=0.01;
2780 Double_t kMinDelayErr=0.0001;
2781
2782 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2783 if (!array) return 0;
2784 TGraphErrors *tlaser=0;
2785 if (array){
2786 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2787 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2788 }
2789 if (!tlaser) return 0;
2790 Int_t npoints0= tlaser->GetN();
2791 if (npoints0==0) return 0;
2792 Double_t *xlaser = new Double_t[npoints0];
2793 Double_t *ylaser = new Double_t[npoints0];
2794 Int_t npoints=0;
2795 for (Int_t i=0;i<npoints0;i++){
2796 //printf("%d\n",i);
2797 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2798 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2799 xlaser[npoints]=tlaser->GetX()[npoints];
2800 ylaser[npoints]=tlaser->GetY()[npoints];
2801 npoints++;
2802 }
2803 //
2804 //
2805 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2806 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2807 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2808 if (index0<0) index0=0;
2809 if (index1>=npoints-1) index1=npoints-1;
f03e487b 2810 if (index1-index0<kMinPoints) {
2811 delete [] ylaser;
2812 delete [] xlaser;
2813 return 0;
2814 }
a23ba1c3 2815 //
2816 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2817 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2818 delete [] ylaser;
2819 delete [] xlaser;
2820 return mean;
2821}
0d1b4cf8 2822
2823
2824
2825
8166d5d7 2826void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
0d1b4cf8 2827 //
2828 // Filter Goofie data
1fabc823 2829 // goofieArray - points will be filtered
2830 // deltaT - smmothing time window
2831 // cutSigma - outler sigma cut in rms
2832 // minVn, maxVd- range absolute cut for variable vd/pt
2833 // - to be tuned
0d1b4cf8 2834 //
0d1b4cf8 2835 // Ignore goofie if not enough points
2836 //
2837 const Int_t kMinPoints = 3;
2838 //
2839
2840 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2841 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2842 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2843 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2844 if (!graphvd) return;
2845 if (graphvd->GetN()<kMinPoints){
2846 delete graphvd;
2847 goofieArray->GetSensorNum(2)->SetGraph(0);
2848 return;
2849 }
2850 //
2851 // 1. Caluclate medians of critical variables
2852 // drift velcocity
2853 // P/T
2854 // area near peak
2855 // area far peak
2856 //
2857 Double_t medianpt=0;
1fabc823 2858 Double_t medianvd=0, sigmavd=0;
0d1b4cf8 2859 Double_t medianan=0;
2860 Double_t medianaf=0;
1fabc823 2861 Int_t entries=graphvd->GetN();
2862 Double_t yvdn[10000];
2863 Int_t nvd=0;
2864 //
2865 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2866 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2867 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2868 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2869 yvdn[nvd++]=graphvd->GetY()[ipoint];
2870 }
2871 if (nvd<kMinPoints){
2872 delete graphvd;
2873 goofieArray->GetSensorNum(2)->SetGraph(0);
2874 return;
2875 }
2876 //
2877 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2878 if (nuni>=kMinPoints){
2879 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2880 }else{
2881 medianvd = TMath::Median(nvd, yvdn);
2882 }
2883
0d1b4cf8 2884 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2885 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2886 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2887 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2888 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2889 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
0d1b4cf8 2890 delete graphpt0;
2891 delete graphpt1;
2892 delete graphan0;
2893 delete graphan1;
2894 delete graphaf0;
2895 delete graphaf1;
2896 //
2897 // 2. Make outlyer graph
2898 //
1fabc823 2899 Int_t nOK=0;
0d1b4cf8 2900 TGraph graphOut(*graphvd);
2901 for (Int_t i=0; i<entries;i++){
2902 //
2903 Bool_t isOut=kFALSE;
1fabc823 2904 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2905 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2906
2907 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2908 isOut|=kTRUE;
0d1b4cf8 2909 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
1fabc823 2910 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2911 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
0d1b4cf8 2912 graphOut.GetY()[i]= (isOut)?1:0;
1fabc823 2913 if (!isOut) nOK++;
0d1b4cf8 2914 }
1fabc823 2915 if (nOK<kMinPoints) {
2916 delete graphvd;
2917 goofieArray->GetSensorNum(2)->SetGraph(0);
2918 return;
2919 }
0d1b4cf8 2920 //
2921 // 3. Filter out outlyers - and smooth
2922 //
2923 TVectorF vmedianArray(goofieArray->NumSensors());
2924 TVectorF vrmsArray(goofieArray->NumSensors());
2925 Double_t xnew[10000];
2926 Double_t ynew[10000];
2927 TObjArray junk;
2928 junk.SetOwner(kTRUE);
2929 Bool_t isOK=kTRUE;
2930 //
2931 //
2932 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2933 isOK=kTRUE;
2934 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2935 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2936 //
2937 if (!sensor) continue;
2938 graphOld = sensor->GetGraph();
2939 if (graphOld) {
2940 sensor->SetGraph(0);
2941 Int_t nused=0;
2942 for (Int_t i=0;i<entries;i++){
2943 if (graphOut.GetY()[i]>0.5) continue;
2944 xnew[nused]=graphOld->GetX()[i];
2945 ynew[nused]=graphOld->GetY()[i];
2946 nused++;
2947 }
2948 graphNew = new TGraph(nused,xnew,ynew);
2949 junk.AddLast(graphNew);
2950 junk.AddLast(graphOld);
2951 Double_t median=0;
2952 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2953 if (graphNew0!=0){
2954 junk.AddLast(graphNew0);
2955 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2956 if (graphNew1!=0){
1fabc823 2957 junk.AddLast(graphNew1);
0d1b4cf8 2958 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2959 if (graphNew2!=0) {
1fabc823 2960 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
0d1b4cf8 2961 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2962 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2963 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
ef1807ef 2964 // AliInfo(Form("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]));
0d1b4cf8 2965 vmedianArray[isensor]=median;
0d1b4cf8 2966 //
2967 }
2968 }
2969 }
2970 }
2971 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2972 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2973 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2974 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2975 (*pcstream)<<"goofieA"<<
2976 Form("isOK_%d.=",isensor)<<isOK<<
2977 Form("s_%d.=",isensor)<<sensor<<
2978 Form("gr_%d.=",isensor)<<graphOld<<
2979 Form("gr0_%d.=",isensor)<<graphNew0<<
2980 Form("gr1_%d.=",isensor)<<graphNew1<<
2981 Form("gr2_%d.=",isensor)<<graphNew2;
1fabc823 2982 if (isOK) sensor->SetGraph(graphNew2);
2983 }
2984 (*pcstream)<<"goofieA"<<
2985 "vmed.="<<&vmedianArray<<
2986 "vrms.="<<&vrmsArray<<
2987 "\n";
2988 junk.Delete(); // delete temoprary graphs
0d1b4cf8 2989
2990}
2991
2992
a980538f 2993
2994
2995
8166d5d7 2996TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
a980538f 2997 //
2998 // Make a statistic matrix
2999 // Input parameters:
3000 // array - TObjArray of AliRelKalmanAlign
3001 // minFraction - minimal ration of accepted tracks
3002 // minStat - minimal statistic (number of accepted tracks)
3003 // maxvd - maximal deviation for the 1
3004 // Output matrix:
3005 // columns - Mean, Median, RMS
3006 // row - parameter type (rotation[3], translation[3], drift[3])
3007 if (!array) return 0;
3008 if (array->GetEntries()<=0) return 0;
cc65e4f5 3009 // Int_t entries = array->GetEntries();
a980538f 3010 Int_t entriesFast = array->GetEntriesFast();
3011 TVectorD state(9);
3012 TVectorD *valArray[9];
3013 for (Int_t i=0; i<9; i++){
3014 valArray[i] = new TVectorD(entriesFast);
3015 }
3016 Int_t naccept=0;
3017 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3018 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3019 if (!kalman) continue;
3020 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3021 if (kalman->GetNUpdates()<minStat) continue;
938ce618 3022 if (Float_t(kalman->GetNUpdates())/Float_t(kalman->GetNTracks())<minFraction) continue;
a980538f 3023 kalman->GetState(state);
3024 for (Int_t ipar=0; ipar<9; ipar++)
3025 (*valArray[ipar])[naccept]=state[ipar];
3026 naccept++;
3027 }
5647625c 3028 //if (naccept<2) return 0;
3029 if (naccept<1) return 0;
a980538f 3030 TMatrixD *pstat=new TMatrixD(9,3);
3031 TMatrixD &stat=*pstat;
3032 for (Int_t ipar=0; ipar<9; ipar++){
3033 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3034 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3035 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3036 }
3037 return pstat;
3038}
3039
3040
8166d5d7 3041TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
a980538f 3042 //
3043 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3044 // Input:
3045 // array - input array
3046 // stat - mean parameters statistic
3047 // direction -
3048 // sigmaCut - maximal allowed deviation from mean in terms of RMS
3049 if (!array) return 0;
3050 if (array->GetEntries()<=0) return 0;
32100b1d 3051 if (!(&stat)) return 0;
cc65e4f5 3052 // error increase in 1 hour
3053 const Double_t kerrsTime[9]={
3054 0.00001, 0.00001, 0.00001,
3055 0.001, 0.001, 0.001,
dfb83639 3056 0.002, 0.01, 0.001};
cc65e4f5 3057 //
3058 //
a980538f 3059 Int_t entries = array->GetEntriesFast();
3060 TObjArray *sArray= new TObjArray(entries);
3061 AliRelAlignerKalman * sKalman =0;
3062 TVectorD state(9);
3063 for (Int_t i=0; i<entries; i++){
3064 Int_t index=(direction)? entries-i-1:i;
3065 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3066 if (!kalman) continue;
3067 Bool_t isOK=kTRUE;
3068 kalman->GetState(state);
3069 for (Int_t ipar=0; ipar<9; ipar++){
3070 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3071 }
3072 if (!sKalman &&isOK) {
3073 sKalman=new AliRelAlignerKalman(*kalman);
3074 sKalman->SetRejectOutliers(kFALSE);
3075 sKalman->SetRunNumber(kalman->GetRunNumber());
3076 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3077 }
3078 if (!sKalman) continue;
3079 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
cc65e4f5 3080 for (Int_t ipar=0; ipar<9; ipar++){
3081// (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3082// (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3083// (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3084 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3085 }
a980538f 3086 sKalman->SetRunNumber(kalman->GetRunNumber());
3087 if (!isOK) sKalman->SetRunNumber(0);
3088 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3089 if (!isOK) continue;
3090 sKalman->SetRejectOutliers(kFALSE);
3091 sKalman->SetRunNumber(kalman->GetRunNumber());
3092 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3093 sKalman->Merge(kalman);
3094 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3095 //sKalman->Print();
3096 }
3097 return sArray;
3098}
3099
8166d5d7 3100TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
cc65e4f5 3101 //
3102 // Merge 2 RelKalman arrays
3103 // Input:
3104 // arrayP - rel kalman in direction plus
3105 // arrayM - rel kalman in direction minus
3106 if (!arrayP) return 0;
3107 if (arrayP->GetEntries()<=0) return 0;
3108 if (!arrayM) return 0;
3109 if (arrayM->GetEntries()<=0) return 0;
3a4b0aed 3110
cc65e4f5 3111 Int_t entries = arrayP->GetEntriesFast();
3112 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3a4b0aed 3113
cc65e4f5 3114 for (Int_t i=0; i<entries; i++){
3115 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3116 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3117 if (!kalmanP) continue;
3118 if (!kalmanM) continue;
3a4b0aed 3119
3120 AliRelAlignerKalman *kalman = NULL;
3121 if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3122 kalman = new AliRelAlignerKalman(*kalmanP);
3123 kalman->Merge(kalmanM);
3124 }
3125 else if (kalmanP->GetRunNumber() == 0) {
3126 kalman = new AliRelAlignerKalman(*kalmanM);
3127 }
3128 else if (kalmanM->GetRunNumber() == 0) {
3129 kalman = new AliRelAlignerKalman(*kalmanP);
3130 }
3131 else
3132 continue;
3133
cc65e4f5 3134 array->AddAt(kalman,i);
3135 }
3a4b0aed 3136
cc65e4f5 3137 return array;
3138}