]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibDButil.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
macro to make the basic calibration trend graphs
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibDButil.cxx
CommitLineData
892226be 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////////
18// //
19// Class providing the calculation of derived quantities (mean,rms,fits,...) //
20// of calibration entries //
21/*
22
23
24*/
25////////////////////////////////////////////////////////////////////////////////
26
27#include <TMath.h>
28#include <TVectorT.h>
29#include <TObjArray.h>
30#include <TGraph.h>
7390f655 31#include <TFile.h>
32#include <TDirectory.h>
949d8707 33#include <TMap.h>
34#include <TGraphErrors.h>
a23ba1c3 35#include <AliCDBStorage.h>
892226be 36#include <AliDCSSensorArray.h>
a23ba1c3 37#include <AliTPCSensorTempArray.h>
892226be 38#include <AliDCSSensor.h>
949d8707 39#include <AliLog.h>
40#include <AliCDBEntry.h>
41#include <AliCDBManager.h>
42#include <AliCDBId.h>
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) ){
93425263 387// printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
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){
463// printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
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
521 printf("Sector %d\n", isec);
522
523 AliTPCCalROC* occupancyROC = fDataQA->GetNoThreshold()->GetCalROC(isec);
524 AliTPCCalROC* nclusterROC = fDataQA->GetNLocalMaxima()->GetCalROC(isec);
525 AliTPCCalROC* qROC = fDataQA->GetMeanCharge()->GetCalROC(isec);
526 AliTPCCalROC* qmaxROC = fDataQA->GetMaxCharge()->GetCalROC(isec);
527 if (!occupancyROC) continue;
528 if (!nclusterROC) continue;
529 if (!qROC) continue;
530 if (!qmaxROC) continue;
531
532 const UInt_t nchannels=occupancyROC->GetNchannels();
533
534 printf("Nchannels %d\n", nchannels);
535
536 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
537
538 vQaOcc[isec] += occupancyROC->GetValue(ichannel);
539 ++normOcc[isec];
540
541 Float_t nClusters = nclusterROC->GetValue(ichannel);
542 normQ[isec] += nClusters;
543 vQaQtot[isec]+=nClusters*qROC->GetValue(ichannel);
544 vQaQmax[isec]+=nClusters*qmaxROC->GetValue(ichannel);
545 }
546 }
547
548 //calculate mean values
549 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
550
551 if (normOcc[isec]>0) vQaOcc[isec] /= normOcc[isec];
552 else vQaOcc[isec] = 0;
553
554 if (normQ[isec]>0) {
555 vQaQtot[isec] /= normQ[isec];
556 vQaQmax[isec] /= normQ[isec];
557 }else {
558
559 vQaQtot[isec] = 0;
560 vQaQmax[isec] = 0;
561 }
562 }
563}
564
892226be 565//_____________________________________________________________________________________
566void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
567{
568 //
569 // Process the Pulser information
570 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
571 //
572
573 const UInt_t infoSize=4;
574 //reset counters to error number
575 vMeanTime.ResizeTo(infoSize);
576 vMeanTime.Zero();
577 //counter
578 TVectorD c(infoSize);
579 //retrieve pulser and ALTRO data
580 if (!fPulserTmean) return;
581 //
582 //get Outliers
583 AliTPCCalROC *rocOut=0x0;
584 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
585 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
586 if (!tmeanROC) continue;
587 rocOut=fPulserOutlier->GetCalROC(isec);
588 UInt_t nchannels=tmeanROC->GetNchannels();
589 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
590 if (rocOut && rocOut->GetValue(ichannel)) continue;
591 Float_t val=tmeanROC->GetValue(ichannel);
592 Int_t type=isec/18;
593 vMeanTime[type]+=val;
594 ++c[type];
595 }
596 }
597 //calculate mean
598 for (UInt_t itype=0; itype<infoSize; ++itype){
599 if (c[itype]>0) vMeanTime[itype]/=c[itype];
600 else vMeanTime[itype]=0;
601 }
602}
603//_____________________________________________________________________________________
604void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
605{
606 //
607 // Get Values from ALTRO configuration data
608 //
609 nMasked=-1;
610 if (!fALTROMasked) return;
611 nMasked=0;
612 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
613 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
614 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
615 if (rocMasked->GetValue(ichannel)) ++nMasked;
616 }
617 }
618}
619//_____________________________________________________________________________________
620void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
621{
622 //
623 // Proces Goofie values, return statistical information of the currently set goofieArray
624 // The meaning of the entries are given below
625 /*
626 1 TPC_ANODE_I_A00_STAT
627 2 TPC_DVM_CO2
628 3 TPC_DVM_DriftVelocity
629 4 TPC_DVM_FCageHV
630 5 TPC_DVM_GainFar
631 6 TPC_DVM_GainNear
632 7 TPC_DVM_N2
633 8 TPC_DVM_NumberOfSparks
634 9 TPC_DVM_PeakAreaFar
635 10 TPC_DVM_PeakAreaNear
636 11 TPC_DVM_PeakPosFar
637 12 TPC_DVM_PeakPosNear
638 13 TPC_DVM_PickupHV
639 14 TPC_DVM_Pressure
640 15 TPC_DVM_T1_Over_P
641 16 TPC_DVM_T2_Over_P
642 17 TPC_DVM_T_Over_P
643 18 TPC_DVM_TemperatureS1
644 */
645 if (!fGoofieArray){
646 Int_t nsensors=19;
647 vecEntries.ResizeTo(nsensors);
648 vecMedian.ResizeTo(nsensors);
649 vecMean.ResizeTo(nsensors);
650 vecRMS.ResizeTo(nsensors);
651 vecEntries.Zero();
652 vecMedian.Zero();
653 vecMean.Zero();
654 vecRMS.Zero();
655 return;
656 }
657 Double_t kEpsilon=0.0000000001;
658 Double_t kBig=100000000000.;
659 Int_t nsensors = fGoofieArray->NumSensors();
660 vecEntries.ResizeTo(nsensors);
661 vecMedian.ResizeTo(nsensors);
662 vecMean.ResizeTo(nsensors);
663 vecRMS.ResizeTo(nsensors);
664 TVectorF values;
665 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
666 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
667 if (gsensor && gsensor->GetGraph()){
668 Int_t npoints = gsensor->GetGraph()->GetN();
669 // filter zeroes
670 values.ResizeTo(npoints);
671 Int_t nused =0;
672 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
673 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
674 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
675 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
676 nused++;
677 }
678 }
679 //
680 vecEntries[isensor]= nused;
681 if (nused>1){
682 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
683 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
684 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
685 }
686 }
687 }
688}
7390f655 689//_____________________________________________________________________________________
690void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
691{
692 //
693 // check the variations of the pedestal data to the reference pedestal data
694 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
695 //
696 const Int_t npar=4;
697 TVectorF vThres(npar); //thresholds
698 Int_t nActive=0; //number of active channels
699
700 //reset and set thresholds
701 pedestalDeviations.ResizeTo(npar);
702 for (Int_t i=0;i<npar;++i){
703 pedestalDeviations.GetMatrixArray()[i]=0;
704 vThres.GetMatrixArray()[i]=(i+1)*.5;
705 }
706 //check all needed data is available
707 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
708 //loop over all channels
709 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
710 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
711 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
712 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
713 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
714 UInt_t nrows=mROC->GetNrows();
715 for (UInt_t irow=0;irow<nrows;++irow){
716 UInt_t npads=mROC->GetNPads(irow);
717 for (UInt_t ipad=0;ipad<npads;++ipad){
718 //don't use masked channels;
719 if (mROC ->GetValue(irow,ipad)) continue;
720 if (mRefROC->GetValue(irow,ipad)) continue;
721 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
722 for (Int_t i=0;i<npar;++i){
723 if (deviation>vThres[i])
724 ++pedestalDeviations.GetMatrixArray()[i];
725 }
726 ++nActive;
727 }//end ipad
728 }//ind irow
729 }//end isec
730 if (nActive>0){
731 for (Int_t i=0;i<npar;++i){
732 pedestalDeviations.GetMatrixArray()[i]/=nActive;
733 }
734 }
735}
736//_____________________________________________________________________________________
737void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
738{
739 //
740 // check the variations of the noise data to the reference noise data
741 // thresholds are 5, 10, 15 and 20 percent respectively.
742 //
743 const Int_t npar=4;
744 TVectorF vThres(npar); //thresholds
745 Int_t nActive=0; //number of active channels
746
747 //reset and set thresholds
748 noiseDeviations.ResizeTo(npar);
749 for (Int_t i=0;i<npar;++i){
750 noiseDeviations.GetMatrixArray()[i]=0;
751 vThres.GetMatrixArray()[i]=(i+1)*.05;
752 }
753 //check all needed data is available
754 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
755 //loop over all channels
756 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
757 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
758 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
759 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
760 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
761 UInt_t nrows=mROC->GetNrows();
762 for (UInt_t irow=0;irow<nrows;++irow){
763 UInt_t npads=mROC->GetNPads(irow);
764 for (UInt_t ipad=0;ipad<npads;++ipad){
765 //don't use masked channels;
766 if (mROC ->GetValue(irow,ipad)) continue;
767 if (mRefROC->GetValue(irow,ipad)) continue;
32100b1d 768 if (nRefROC->GetValue(irow,ipad)==0) continue;
7390f655 769 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
770 for (Int_t i=0;i<npar;++i){
771 if (deviation>vThres[i])
772 ++noiseDeviations.GetMatrixArray()[i];
773 }
774 ++nActive;
775 }//end ipad
776 }//ind irow
777 }//end isec
778 if (nActive>0){
779 for (Int_t i=0;i<npar;++i){
780 noiseDeviations.GetMatrixArray()[i]/=nActive;
781 }
782 }
783}
784//_____________________________________________________________________________________
785void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
786 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
787{
788 //
789 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
790 // thresholds are .5, 1, 5 and 10 percent respectively.
791 //
792 //
793 const Int_t npar=4;
794 TVectorF vThres(npar); //thresholds
795 Int_t nActive=0; //number of active channels
796
797 //reset and set thresholds
798 pulserQdeviations.ResizeTo(npar);
799 for (Int_t i=0;i<npar;++i){
800 pulserQdeviations.GetMatrixArray()[i]=0;
801 }
802 npadsOutOneTB=0;
803 npadsOffAdd=0;
804 varQMean=0;
805 vThres.GetMatrixArray()[0]=.005;
806 vThres.GetMatrixArray()[1]=.01;
807 vThres.GetMatrixArray()[2]=.05;
808 vThres.GetMatrixArray()[3]=.1;
809 //check all needed data is available
810 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
811 //
812 UpdateRefPulserOutlierMap();
813 //loop over all channels
732e90a8 814 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 815 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
816 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
817 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
732e90a8 818// AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
7390f655 819 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
820 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
821 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
d30aa177 822 Float_t ptmean=ptROC->GetMean(oROC);
7390f655 823 UInt_t nrows=mROC->GetNrows();
824 for (UInt_t irow=0;irow<nrows;++irow){
825 UInt_t npads=mROC->GetNPads(irow);
826 for (UInt_t ipad=0;ipad<npads;++ipad){
827 //don't use masked channels;
828 if (mROC ->GetValue(irow,ipad)) continue;
829 if (mRefROC->GetValue(irow,ipad)) continue;
830 //don't user edge pads
831 if (ipad==0||ipad==npads-1) continue;
832 //data
833 Float_t pq=pqROC->GetValue(irow,ipad);
834 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
835 Float_t pt=ptROC->GetValue(irow,ipad);
836// Float_t ptRef=ptRefROC->GetValue(irow,ipad);
837 //comparisons q
838 Float_t deviation=TMath::Abs(pq/pqRef-1);
839 for (Int_t i=0;i<npar;++i){
840 if (deviation>vThres[i])
841 ++pulserQdeviations.GetMatrixArray()[i];
842 }
843 if (pqRef>11&&pq<11) ++npadsOffAdd;
844 varQMean+=pq-pqRef;
845 //comparisons t
d30aa177 846 if (TMath::Abs(pt-ptmean)>1) ++npadsOutOneTB;
7390f655 847 ++nActive;
848 }//end ipad
849 }//ind irow
850 }//end isec
851 if (nActive>0){
852 for (Int_t i=0;i<npar;++i){
853 pulserQdeviations.GetMatrixArray()[i]/=nActive;
854 varQMean/=nActive;
855 }
856 }
857}
892226be 858//_____________________________________________________________________________________
859void AliTPCcalibDButil::UpdatePulserOutlierMap()
7390f655 860{
861 //
8166d5d7 862 // Update the outlier map of the pulser data
7390f655 863 //
864 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
865}
866//_____________________________________________________________________________________
867void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
868{
869 //
8166d5d7 870 // Update the outlier map of the pulser reference data
7390f655 871 //
872 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
873}
874//_____________________________________________________________________________________
875void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
892226be 876{
877 //
878 // Create a map that contains outliers from the Pulser calibration data.
879 // The outliers include masked channels, edge pads and pads with
880 // too large timing and charge variations.
7390f655 881 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
892226be 882 // those do not contain masked and edge pads
883 //
7390f655 884 if (!pulT||!pulQ) {
892226be 885 //reset map
7390f655 886 pulOut->Multiply(0.);
892226be 887 fNpulserOutliers=-1;
888 return;
889 }
890 AliTPCCalROC *rocMasked=0x0;
891 fNpulserOutliers=0;
892
893 //Create Outlier Map
894 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 895 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
896 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
897 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
892226be 898 if (!tmeanROC||!qmeanROC) {
899 //reset outliers in this ROC
900 outROC->Multiply(0.);
901 continue;
902 }
903 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
904// Double_t dummy=0;
905// Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
906// Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
907 UInt_t nrows=tmeanROC->GetNrows();
908 for (UInt_t irow=0;irow<nrows;++irow){
909 UInt_t npads=tmeanROC->GetNPads(irow);
910 for (UInt_t ipad=0;ipad<npads;++ipad){
911 Int_t outlier=0,masked=0;
912 Float_t q=qmeanROC->GetValue(irow,ipad);
913 Float_t t=tmeanROC->GetValue(irow,ipad);
914 //masked channels are outliers
915 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
916 //edge pads are outliers
917 if (ipad==0||ipad==npads-1) masked=1;
918 //channels with too large charge or timing deviation from the meadian are outliers
919// if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
920 if (q<fPulQminLimit && !masked) outlier=1;
921 //check for nan
922 if ( !(q<10000000) || !(t<10000000)) outlier=1;
923 outROC->SetValue(irow,ipad,outlier+masked);
924 fNpulserOutliers+=outlier;
925 }
926 }
927 }
928}
929//_____________________________________________________________________________________
2cb269df 930AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
892226be 931{
932 //
933 // Create pad time0 object from pulser and/or CE data, depending on the selected model
934 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
935 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
732e90a8 936 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
892226be 937 //
2cb269df 938 // In case model 2 is invoked - gy arival time gradient is also returned
939 //
940 gyA=0;
941 gyC=0;
892226be 942 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
943 // decide between different models
944 if (model==0||model==1){
945 TVectorD vMean;
946 if (model==1) ProcessPulser(vMean);
947 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
948 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
949 if (!rocPulTmean) continue;
950 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
951 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
952 Float_t mean=rocPulTmean->GetMean(rocOut);
953 //treat case where a whole partition is masked
8166d5d7 954 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
892226be 955 if (model==1) {
956 Int_t type=isec/18;
957 mean=vMean[type];
958 }
959 UInt_t nrows=rocTime0->GetNrows();
960 for (UInt_t irow=0;irow<nrows;++irow){
961 UInt_t npads=rocTime0->GetNPads(irow);
962 for (UInt_t ipad=0;ipad<npads;++ipad){
963 Float_t time=rocPulTmean->GetValue(irow,ipad);
964 //in case of an outlier pad use the mean of the altro values.
965 //This should be the most precise guess in that case.
966 if (rocOut->GetValue(irow,ipad)) {
967 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
8166d5d7 968 if ( TMath::Abs(time)<kAlmost0 ) time=mean;
892226be 969 }
970 Float_t val=time-mean;
971 rocTime0->SetValue(irow,ipad,val);
972 }
973 }
974 }
2cb269df 975 } else if (model==2){
976 Double_t pgya,pgyc,pchi2a,pchi2c;
977 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
978 fCETmean->Add(padPulser,-1.);
732e90a8 979 TVectorD vA,vC;
980 AliTPCCalPad outCE("outCE","outCE");
981 Int_t nOut;
2cb269df 982 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
983 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
732e90a8 984// AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
2cb269df 985 if (!padFit) { delete padPulser; return 0;}
986 gyA=vA[2];
987 gyC=vC[2];
988 fCETmean->Add(padPulser,1.);
732e90a8 989 padTime0->Add(fCETmean);
2cb269df 990 padTime0->Add(padFit,-1);
991 delete padPulser;
732e90a8 992 TVectorD vFitROC;
993 TMatrixD mFitROC;
994 Float_t chi2;
995 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
996 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
997 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
998 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
999 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
1000 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
1001 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
1002 Float_t mean=rocPulTmean->GetMean(rocOutPul);
8166d5d7 1003 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
732e90a8 1004 UInt_t nrows=rocTime0->GetNrows();
1005 for (UInt_t irow=0;irow<nrows;++irow){
1006 UInt_t npads=rocTime0->GetNPads(irow);
1007 for (UInt_t ipad=0;ipad<npads;++ipad){
1008 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
1009 if (rocOutCE->GetValue(irow,ipad)){
1010 Float_t valOut=rocCEfit->GetValue(irow,ipad);
1011 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
1012 rocTime0->SetValue(irow,ipad,valOut);
1013 }
1014 }
1015 }
1016 delete rocCEfit;
1017 }
1018 delete padFit;
892226be 1019 }
2cb269df 1020 Double_t median = padTime0->GetMedian();
1021 padTime0->Add(-median); // normalize to median
892226be 1022 return padTime0;
1023}
1024//_____________________________________________________________________________________
98a4cc77 1025Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *const rocOut)
892226be 1026{
8166d5d7 1027 //
1028 // GetMeanAlto information
1029 //
892226be 1030 if (roc==0) return 0.;
1031 const Int_t sector=roc->GetSector();
1032 AliTPCROC *tpcRoc=AliTPCROC::Instance();
1033 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
1034 Float_t mean=0;
1035 Int_t n=0;
1036
1037 //loop over a small range around the requested pad (+-10 rows/pads)
1038 for (Int_t irow=row-10;irow<row+10;++irow){
1039 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
1040 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
1041 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
1042 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
1043 if (altroRoc!=altroCurr) continue;
1044 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
1045 Float_t val=roc->GetValue(irow,ipad);
1046 mean+=val;
1047 ++n;
1048 }
1049 }
1050 if (n>0) mean/=n;
1051 return mean;
1052}
7390f655 1053//_____________________________________________________________________________________
1054void AliTPCcalibDButil::SetRefFile(const char* filename)
1055{
1056 //
1057 // load cal pad objects form the reference file
1058 //
1059 TDirectory *currDir=gDirectory;
1060 TFile f(filename);
1061 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
1062 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
1063 //pulser data
1064 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
1065 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
1066 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
1067 //CE data
1068 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
1069 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
1070 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
1071 //Altro data
1072// fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
1073// fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
1074// fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
1075// fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
1076 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
1077 f.Close();
1078 currDir->cd();
1079}
949d8707 1080//_____________________________________________________________________________________
1081void AliTPCcalibDButil::UpdateRefDataFromOCDB()
1082{
1083 //
1084 // set reference data from OCDB Reference map
1085 //
1086 if (!fRefMap) {
1087 AliWarning("Referenc map not set!");
1088 return;
1089 }
1090
1091 TString cdbPath;
1092 AliCDBEntry* entry = 0x0;
1093 Bool_t hasAnyChanged=kFALSE;
1094
1095 //pedestals
1096 cdbPath="TPC/Calib/Pedestals";
1097 if (HasRefChanged(cdbPath.Data())){
1098 hasAnyChanged=kTRUE;
1099 //delete old entries
1100 if (fRefPedestals) delete fRefPedestals;
1101 if (fRefPedestalMasked) delete fRefPedestalMasked;
1102 fRefPedestals=fRefPedestalMasked=0x0;
1103 //get new entries
1104 entry=GetRefEntry(cdbPath.Data());
1105 if (entry){
1106 entry->SetOwner(kTRUE);
1107 fRefPedestals=GetRefCalPad(entry);
1108 delete entry;
1109 fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1110 }
1111 }
7390f655 1112
949d8707 1113 //noise
1114 cdbPath="TPC/Calib/PadNoise";
1115 if (HasRefChanged(cdbPath.Data())){
1116 hasAnyChanged=kTRUE;
1117 //delete old entry
1118 if (fRefPadNoise) delete fRefPadNoise;
1119 fRefPadNoise=0x0;
1120 //get new entry
1121 entry=GetRefEntry(cdbPath.Data());
1122 if (entry){
1123 entry->SetOwner(kTRUE);
1124 fRefPadNoise=GetRefCalPad(entry);
1125 delete entry;
1126 }
1127 }
1128
1129 //pulser
1130 cdbPath="TPC/Calib/Pulser";
1131 if (HasRefChanged(cdbPath.Data())){
1132 hasAnyChanged=kTRUE;
1133 //delete old entries
1134 if (fRefPulserTmean) delete fRefPulserTmean;
1135 if (fRefPulserTrms) delete fRefPulserTrms;
1136 if (fRefPulserQmean) delete fRefPulserQmean;
1137 if (fRefPulserMasked) delete fRefPulserMasked;
1138 fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1139 //get new entries
1140 entry=GetRefEntry(cdbPath.Data());
1141 if (entry){
1142 entry->SetOwner(kTRUE);
1143 fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1144 fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1145 fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1146 delete entry;
1147 fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1148 }
1149 }
2cb269df 1150
949d8707 1151 //ce
1152 cdbPath="TPC/Calib/CE";
1153 if (HasRefChanged(cdbPath.Data())){
1154 hasAnyChanged=kTRUE;
1155 //delete old entries
1156 if (fRefCETmean) delete fRefCETmean;
1157 if (fRefCETrms) delete fRefCETrms;
1158 if (fRefCEQmean) delete fRefCEQmean;
1159 if (fRefCEMasked) delete fRefCEMasked;
1160 fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1161 //get new entries
1162 entry=GetRefEntry(cdbPath.Data());
1163 if (entry){
1164 entry->SetOwner(kTRUE);
1165 fRefCETmean=GetRefCalPad(entry,"CETmean");
1166 fRefCETrms=GetRefCalPad(entry,"CETrms");
1167 fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1168 delete entry;
1169 fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1170 }
1171 }
1172
1173 //altro data
1174 cdbPath="TPC/Calib/AltroConfig";
1175 if (HasRefChanged(cdbPath.Data())){
1176 hasAnyChanged=kTRUE;
1177 //delete old entries
1178 if (fRefALTROFPED) delete fRefALTROFPED;
1179 if (fRefALTROZsThr) delete fRefALTROZsThr;
1180 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1181 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1182 if (fRefALTROMasked) delete fRefALTROMasked;
1183 fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1184 //get new entries
1185 entry=GetRefEntry(cdbPath.Data());
1186 if (entry){
1187 entry->SetOwner(kTRUE);
1188 fRefALTROFPED=GetRefCalPad(entry,"FPED");
1189 fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1190 fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1191 fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1192 fRefALTROMasked=GetRefCalPad(entry,"Masked");
1193 delete entry;
1194 }
1195 }
1196
1197 //raw data
1198 /*
1199 cdbPath="TPC/Calib/Raw";
1200 if (HasRefChanged(cdbPath.Data())){
1201 hasAnyChanged=kTRUE;
1202 //delete old entry
1203 if (fRefCalibRaw) delete fRefCalibRaw;
1204 //get new entry
1205 entry=GetRefEntry(cdbPath.Data());
1206 if (entry){
1207 entry->SetOwner(kTRUE);
1208 TObjArray *arr=(TObjArray*)entry->GetObject();
1209 if (!arr){
a980538f 1210 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1211 } else {
1212 fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1213 }
1214 }
1215 }
1216 */
2cb269df 1217
949d8707 1218 //data qa
1219 cdbPath="TPC/Calib/QA";
1220 if (HasRefChanged(cdbPath.Data())){
1221 hasAnyChanged=kTRUE;
1222 //delete old entry
1223 if (fRefDataQA) delete fRefDataQA;
1224 //get new entry
1225 entry=GetRefEntry(cdbPath.Data());
1226 if (entry){
1227 entry->SetOwner(kTRUE);
abb20887 1228 fRefDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1229 if (!fRefDataQA){
a980538f 1230 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1231 } else {
668e6ebb 1232 fRefDataQA=(AliTPCdataQA*)fRefDataQA->Clone();
949d8707 1233 }
1234 delete entry;
1235 }
1236 }
1237
1238
1239//update current reference maps
1240 if (hasAnyChanged){
1241 if (fCurrentRefMap) delete fCurrentRefMap;
1242 fCurrentRefMap=(TMap*)fRefMap->Clone();
1243 }
1244}
1245//_____________________________________________________________________________________
1246AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1247{
1248 //
1249 // TObjArray object type case
1250 // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1251 //
1252 AliTPCCalPad *pad=0x0;
1253 TObjArray *arr=(TObjArray*)entry->GetObject();
1254 if (!arr){
a980538f 1255 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1256 return pad;
1257 }
1258 pad=(AliTPCCalPad*)arr->FindObject(objName);
1259 if (!pad) {
a980538f 1260 AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
949d8707 1261 return pad;
1262 }
1263 return (AliTPCCalPad*)pad->Clone();
1264}
1265//_____________________________________________________________________________________
1266AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1267{
1268 //
1269 // AliTPCCalPad object type case
1270 // cast object to a calPad and store it in 'pad'
1271 //
1272 AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1273 if (!pad) {
a980538f 1274 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1275 return 0x0;
1276 }
1277 pad=(AliTPCCalPad*)pad->Clone();
1278 return pad;
1279}
1280//_____________________________________________________________________________________
1281AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1282{
1283 //
1284 // set altro masked channel map for 'cdbPath'
1285 //
1286 AliTPCCalPad* pad=0x0;
1287 const Int_t run=GetReferenceRun(cdbPath);
1288 if (run<0) {
a980538f 1289 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1290 return pad;
1291 }
1292 AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1293 if (!entry) {
a980538f 1294 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1295 return pad;
1296 }
1297 pad=GetRefCalPad(entry,"Masked");
1298 if (pad) pad->SetNameTitle(name,name);
1299 entry->SetOwner(kTRUE);
1300 delete entry;
1301 return pad;
1302}
1303//_____________________________________________________________________________________
1304void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1305 //
1306 // Get Reference map
1307 //
1308 if (run<0) run=fCalibDB->GetRun();
1309 TString cdbPath="TPC/Calib/Ref";
1310 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1311 if (!entry) {
a980538f 1312 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
949d8707 1313 fRefMap=0;
1314 return;
1315 }
1316 entry->SetOwner(kTRUE);
1317 fRefMap=(TMap*)(entry->GetObject());
1318 AliCDBId &id=entry->GetId();
1319 fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1320}
1321//_____________________________________________________________________________________
1322Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1323{
1324 //
1325 // check whether a reference cdb entry has changed
1326 //
1327 if (!fCurrentRefMap) return kTRUE;
1328 if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1329 return kFALSE;
1330}
1331//_____________________________________________________________________________________
1332AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1333{
1334 //
1335 // get the reference AliCDBEntry for 'cdbPath'
1336 //
1337 const Int_t run=GetReferenceRun(cdbPath);
1338 if (run<0) {
a980538f 1339 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1340 return 0;
1341 }
1342 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1343 if (!entry) {
a980538f 1344 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1345 return 0;
1346 }
1347 return entry;
1348}
1349//_____________________________________________________________________________________
26d8859c 1350Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
949d8707 1351 //
1352 // Get reference run number for the specified OCDB path
1353 //
1354 if (!fCurrentRefMap) return -2;
1355 TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1356 if (!str) return -2;
58affc42 1357 return (Int_t)str->GetString().Atoi();
949d8707 1358}
1359//_____________________________________________________________________________________
26d8859c 1360Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
949d8707 1361 //
1362 // Get reference run number for the specified OCDB path
1363 //
1364 if (!fRefMap) return -1;
1365 TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1366 if (!str) return -1;
58affc42 1367 return (Int_t)str->GetString().Atoi();
949d8707 1368}
1369//_____________________________________________________________________________________
8166d5d7 1370AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
2cb269df 1371 //
1372 // Author: marian.ivanov@cern.ch
1373 //
1374 // Create outlier map for CE study
1375 // Parameters:
1376 // Return value - outlyer map
1377 // noutlyersCE - number of outlyers
1378 // minSignal - minimal total Q signal
1379 // cutRMSMin - minimal width of the signal in respect to the median
1380 // cutRMSMax - maximal width of the signal in respect to the median
1381 // cutMaxDistT - maximal deviation from time median per chamber
1382 //
1383 // Outlyers criteria:
1384 // 0. Exclude masked pads
1385 // 1. Exclude first two rows in IROC and last two rows in OROC
1386 // 2. Exclude edge pads
1387 // 3. Exclude channels with too large variations
1388 // 4. Exclude pads with too small signal
1389 // 5. Exclude signal with outlyers RMS
1390 // 6. Exclude channels to far from the chamber median
1391 noutliersCE=0;
1392 //create outlier map
1393 AliTPCCalPad *out=ceOut;
1394 if (!out) out= new AliTPCCalPad("outCE","outCE");
1395 AliTPCCalROC *rocMasked=0x0;
1396 if (!fCETmean) return 0;
1397 if (!fCETrms) return 0;
1398 if (!fCEQmean) return 0;
1399 //
1400 //loop over all channels
1401 //
1402 Double_t rmsMedian = fCETrms->GetMedian();
1403 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1404 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
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 //
f03e487b 1411 if (!rocData || !rocCEQ || !rocCETrms) {
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 }
1779 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
aec07590 1780 if (nused<kMinPoints) {
f03e487b 1781 delete [] tdelta;
aec07590 1782 printf("AliFatal: No time offset calibration available\n");
1783 return 0;
1784 }
817766d5 1785 Double_t median = TMath::Median(nused,tdelta);
1786 Double_t mean = TMath::Mean(nused,tdelta);
a9f487a1 1787 delete [] tdelta;
817766d5 1788 return (valType==0) ? median:mean;
1789}
1790
a23ba1c3 1791Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
817766d5 1792 //
1793 // Get the correction of the drift velocity
1794 // combining information from the laser track calibration
1795 // and from cosmic calibration
1796 //
a23ba1c3 1797 // dist - return value - distance to closest point in graph
817766d5 1798 // run - run number
1799 // timeStamp - tim stamp in seconds
1800 // deltaT - integration period to calculate time0 offset
1801 // deltaTLaser -max validity of laser data
1802 // valType - 0 - median, 1- mean
1803 //
1804 // Integration vaues are just recomendation - if not possible to get points
1805 // automatically increase the validity by factor 2
1806 // (recursive algorithm until one month of data taking)
1807 //
1808 //
1809 //
1810 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1811 if (!array) {
1812 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1813 }
1814 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1815 if (!array) return 0;
1816 TGraphErrors *cosmicAll=0;
1817 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1818 if (!cosmicAll) return 0;
a23ba1c3 1819 Double_t grY=0;
1820 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1821
817766d5 1822 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
5647625c 1823 Double_t vcosmic = AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1e722a63 1824 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
817766d5 1825 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
cc65e4f5 1826 return vcosmic-t0;
817766d5 1827
1828 /*
1829 Example usage:
1830
1831 Int_t run=89000
1832 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1833 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1834 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1835 //
1836 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1837 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1838 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1839 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1840
1841 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1842 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1843
1844 */
1845
1846}
1847
949d8707 1848const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1849{
1850 //
1851 // Create a default name for the gui file
1852 //
1853
1854 return Form("guiRefTreeRun%s.root",GetRefValidity());
1855}
1856
1857Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1858{
1859 //
1860 // Create a gui reference tree
1861 // if dirname and filename are empty default values will be used
1862 // this is the recommended way of using this function
1863 // it allows to check whether a file with the given run validity alredy exists
1864 //
1865 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1866 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1867 return kFALSE;
1868 }
1869
1870 TString file=filename;
1871 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1872
1873 AliTPCPreprocessorOnline prep;
1874 //noise and pedestals
1875 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1876 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1877 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1878 //pulser data
1879 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1880 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1881 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1882 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1883 //CE data
1884 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1885 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1886 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1887 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1888 //Altro data
1889 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1890 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1891 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1892 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1893 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1894 //QA
1895 AliTPCdataQA *dataQA=fRefDataQA;
1896 if (dataQA) {
1897 if (dataQA->GetNLocalMaxima())
1898 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1899 if (dataQA->GetMaxCharge())
1900 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1901 if (dataQA->GetMeanCharge())
1902 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1903 if (dataQA->GetNoThreshold())
1904 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1905 if (dataQA->GetNTimeBins())
1906 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1907 if (dataQA->GetNPads())
1908 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1909 if (dataQA->GetTimePosition())
1910 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1911 }
1912 prep.DumpToFile(file.Data());
1913 return kTRUE;
1914}
1915
a23ba1c3 1916Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1917 //
78f17711 1918 // Get the correction of the drift velocity using the offline laser tracks calbration
a23ba1c3 1919 //
1920 // run - run number
1921 // timeStamp - tim stamp in seconds
1922 // deltaT - integration period to calculate time0 offset
1923 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1924 // Note in case no data form both A and C side - the value from active side used
1925 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
78f17711 1926
1927 return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1928}
1929
1930Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1931 //
1932 // Get the correction of the drift velocity using the online laser tracks calbration
1933 //
1934 // run - run number
1935 // timeStamp - tim stamp in seconds
1936 // deltaT - integration period to calculate time0 offset
1937 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1938 // Note in case no data form both A and C side - the value from active side used
1939 TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1940
1941 Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1942 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1943 if (!param) return 0;
1944
1945 //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1946 dv*=param->GetDriftV()/2.61301900000000000e+06;
1947 if (dv>1e-20) dv=1/dv-1;
1948 else return 0;
1949 // T/P correction
1950 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1951
1952 AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
2cd3a697 1953 AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
78f17711 1954
1955 Double_t corrPTA=0;
1956 Double_t corrPTC=0;
1957
1958 if (temp&&press) {
1959 AliTPCCalibVdrift corr(temp,press,0);
1960 corrPTA=corr.GetPTRelative(timeStamp,0);
1961 corrPTC=corr.GetPTRelative(timeStamp,1);
1962 }
1963
1964 if (side==0) dv -= corrPTA;
1965 if (side==1) dv -= corrPTC;
1966 if (side==2) dv -= (corrPTA+corrPTC)/2;
1967
1968 return dv;
1969}
1970
1971Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1972 Int_t side, TObjArray * const array){
1973 //
1974 // common drift velocity retrieval for online and offline method
1975 //
a23ba1c3 1976 TGraphErrors *grlaserA=0;
1977 TGraphErrors *grlaserC=0;
1978 Double_t vlaserA=0, vlaserC=0;
1979 if (!array) return 0;
1980 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1981 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1982 Double_t deltaY;
5647625c 1983 if (grlaserA && grlaserA->GetN()>0) {
a23ba1c3 1984 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1985 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1986 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1987 }
5647625c 1988 if (grlaserC && grlaserC->GetN()>0) {
a23ba1c3 1989 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1990 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1991 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1992 }
1993 if (side==0) return vlaserA;
1994 if (side==1) return vlaserC;
1995 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1996 if (!grlaserA) return vlaserC;
1997 if (!grlaserC) return vlaserA;
1998 return mdrift;
1999}
2000
2001
2002Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
2003 //
2004 // Get the correction of the drift velocity using the CE laser data
2005 // combining information from the CE, laser track calibration
2006 // and P/T calibration
2007 //
2008 // run - run number
2009 // timeStamp - tim stamp in seconds
2010 // deltaT - integration period to calculate time0 offset
2011 // side - 0 - A side, 1 - C side, 2 - mean from both sides
2012 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
0d1b4cf8 2013 if (!arrT) return 0;
a23ba1c3 2014 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
2015 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2016 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2017 //
2018 //
2019 Double_t corrPTA = 0, corrPTC=0;
2020 Double_t ltime0A = 0, ltime0C=0;
2021 Double_t gry=0;
2022 Double_t corrA=0, corrC=0;
2023 Double_t timeA=0, timeC=0;
5647625c 2024 const Double_t kEpsilon = 0.00001;
a23ba1c3 2025 TGraph *graphA = (TGraph*)arrT->At(72);
2026 TGraph *graphC = (TGraph*)arrT->At(73);
2027 if (!graphA && !graphC) return 0.;
2028 if (graphA &&graphA->GetN()>0) {
2029 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2030 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2031 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2032 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2033 if(ltime0A < kEpsilon) return 0;
a23ba1c3 2034 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2035 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2036 corrA-=corrPTA;
2037 }
2038 if (graphC&&graphC->GetN()>0){
2039 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2040 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2041 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2042 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2043 if(ltime0C < kEpsilon) return 0;
2044if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2045 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2046 corrC-=corrPTC;
2047 }
2048
2049 if (side ==0 ) return corrA;
2050 if (side ==1 ) return corrC;
2051 Double_t corrM= (corrA+corrC)*0.5;
2052 if (!graphA) corrM=corrC;
2053 if (!graphC) corrM=corrA;
2054 return corrM;
2055}
2056
cc65e4f5 2057Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2058 //
2059 // return drift velocity using the TPC-ITS matchin method
2060 // return also distance to the closest point
2061 //
2062 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2063 TGraphErrors *graph=0;
2064 dist=0;
2065 if (!array) return 0;
5647625c 2066 //array->ls();
cc65e4f5 2067 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
2068 if (!graph) return 0;
2069 Double_t deltaY;
2070 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2071 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2072 return value;
2073}
2074
2075Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2076 //
2077 // Get time dependent time 0 (trigger delay in cm) correction
2078 // Arguments:
2079 // timestamp - timestamp
2080 // run - run number
2081 //
2082 // Notice - Extrapolation outside of calibration range - using constant function
2083 //
2084 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2085 TGraphErrors *graph=0;
2086 dist=0;
2087 if (!array) return 0;
2088 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
2089 if (!graph) return 0;
2090 Double_t deltaY;
2091 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2092 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2093 return value;
2094}
2095
2096
a23ba1c3 2097
2098
2099
2100Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2101 //
2102 // VERY obscure method - we need something in framework
2103 // Find the TPC runs with temperature OCDB entry
2104 // cache the start and end of the run
2105 //
2106 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2107 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2108 if (!storage) return 0;
2109 TString path=storage->GetURI();
2110 TString runsT;
2111 {
2112 TString command;
2113 if (path.Contains("local")){ // find the list if local system
2114 path.ReplaceAll("local://","");
2115 path+="TPC/Calib/Temperature";
2116 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2117 }
2118 runsT=gSystem->GetFromPipe(command);
2119 }
2120 TObjArray *arr= runsT.Tokenize("r");
2121 if (!arr) return 0;
2122 //
2123 TArrayI indexes(arr->GetEntries());
2124 TArrayI runs(arr->GetEntries());
2125 Int_t naccept=0;
2126 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2127 Int_t irunN = atoi(arr->At(irun)->GetName());
2128 if (irunN<startRun) continue;
2129 if (irunN>stopRun) continue;
2130 runs[naccept]=irunN;
2131 naccept++;
2132 }}
2133 fRuns.Set(naccept);
2134 fRunsStart.Set(fRuns.fN);
2135 fRunsStop.Set(fRuns.fN);
2136 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2137 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2138
2139 //
2140 AliCDBEntry * entry = 0;
2141 {for (Int_t irun=0;irun<fRuns.fN; irun++){
2142 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2143 if (!entry) continue;
2144 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2145 if (!tmpRun) continue;
2146 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2147 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2148 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
2149 }}
2150 return fRuns.fN;
2151}
2152
2153
2154Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2155 //
2156 // binary search - find the run for given time stamp
2157 //
2158 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2159 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2160 Int_t cindex = -1;
2161 for (Int_t index=index0; index<=index1; index++){
2162 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2163 if (debug) {
2164 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2165 }
2166 }
2167 if (cindex<0) cindex =(index0+index1)/2;
2168 if (cindex<0) {
2169 return 0;
2170 }
2171 return fRuns[cindex];
2172}
2173
2174
2175
2176
2177
2178TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2179 //
2180 // filter outlyer measurement
2181 // Only points around median +- sigmaCut filtered
2182 if (!graph) return 0;
2183 Int_t kMinPoints=2;
2184 Int_t npoints0 = graph->GetN();
2185 Int_t npoints=0;
2186 Float_t rmsY=0;
a23ba1c3 2187 //
2188 //
2189 if (npoints0<kMinPoints) return 0;
a9f487a1 2190
2191 Double_t *outx=new Double_t[npoints0];
2192 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2193 for (Int_t iter=0; iter<3; iter++){
2194 npoints=0;
2195 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2196 if (graph->GetY()[ipoint]==0) continue;
2197 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2198 outx[npoints] = graph->GetX()[ipoint];
2199 outy[npoints] = graph->GetY()[ipoint];
2200 npoints++;
2201 }
2202 if (npoints<=1) break;
2203 medianY =TMath::Median(npoints,outy);
2204 rmsY =TMath::RMS(npoints,outy);
2205 }
2206 TGraph *graphOut=0;
2207 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2208 delete [] outx;
2209 delete [] outy;
a23ba1c3 2210 return graphOut;
2211}
2212
2213
2214TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2215 //
2216 // filter outlyer measurement
2217 // Only points around median +- cut filtered
2218 if (!graph) return 0;
2219 Int_t kMinPoints=2;
2220 Int_t npoints0 = graph->GetN();
2221 Int_t npoints=0;
2222 Float_t rmsY=0;
a23ba1c3 2223 //
2224 //
2225 if (npoints0<kMinPoints) return 0;
a9f487a1 2226
2227 Double_t *outx=new Double_t[npoints0];
2228 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2229 for (Int_t iter=0; iter<3; iter++){
2230 npoints=0;
2231 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2232 if (graph->GetY()[ipoint]==0) continue;
2233 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2234 outx[npoints] = graph->GetX()[ipoint];
2235 outy[npoints] = graph->GetY()[ipoint];
2236 npoints++;
2237 }
2238 if (npoints<=1) break;
2239 medianY =TMath::Median(npoints,outy);
2240 rmsY =TMath::RMS(npoints,outy);
2241 }
2242 TGraph *graphOut=0;
2243 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2244 delete [] outx;
2245 delete [] outy;
a23ba1c3 2246 return graphOut;
2247}
2248
2249
2250
8166d5d7 2251TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
a23ba1c3 2252 //
2253 // filter outlyer measurement
2254 // Only points with normalized errors median +- sigmaCut filtered
2255 //
2256 Int_t kMinPoints=10;
2257 Int_t npoints0 = graph->GetN();
2258 Int_t npoints=0;
2259 Float_t medianErr=0, rmsErr=0;
a9f487a1 2260 //
2261 //
2262 if (npoints0<kMinPoints) return 0;
2263
a23ba1c3 2264 Double_t *outx=new Double_t[npoints0];
2265 Double_t *outy=new Double_t[npoints0];
2266 Double_t *erry=new Double_t[npoints0];
2267 Double_t *nerry=new Double_t[npoints0];
2268 Double_t *errx=new Double_t[npoints0];
a9f487a1 2269
a23ba1c3 2270 for (Int_t iter=0; iter<3; iter++){
2271 npoints=0;
2272 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2273 nerry[npoints] = graph->GetErrorY(ipoint);
2274 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2275 erry[npoints] = graph->GetErrorY(ipoint);
2276 outx[npoints] = graph->GetX()[ipoint];
2277 outy[npoints] = graph->GetY()[ipoint];
2278 errx[npoints] = graph->GetErrorY(ipoint);
2279 npoints++;
2280 }
2281 if (npoints==0) break;
2282 medianErr=TMath::Median(npoints,erry);
2283 medianY =TMath::Median(npoints,outy);
2284 rmsErr =TMath::RMS(npoints,erry);
2285 }
2286 TGraphErrors *graphOut=0;
2287 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2288 delete []outx;
2289 delete []outy;
a23ba1c3 2290 delete []erry;
a9f487a1 2291 delete []nerry;
2292 delete []errx;
a23ba1c3 2293 return graphOut;
2294}
2295
2296
2297void AliTPCcalibDButil::Sort(TGraph *graph){
2298 //
2299 // sort array - neccessay for approx
2300 //
2301 Int_t npoints = graph->GetN();
2302 Int_t *indexes=new Int_t[npoints];
2303 Double_t *outx=new Double_t[npoints];
2304 Double_t *outy=new Double_t[npoints];
2305 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2306 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2307 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2308 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2309 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2310
a9f487a1 2311 delete [] indexes;
2312 delete [] outx;
2313 delete [] outy;
a23ba1c3 2314}
2315void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2316 //
2317 // smmoth graph - mean on the interval
2318 //
2319 Sort(graph);
2320 Int_t npoints = graph->GetN();
2321 Double_t *outy=new Double_t[npoints];
0d1b4cf8 2322
a23ba1c3 2323 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2324 Double_t lx=graph->GetX()[ipoint];
2325 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2326 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2327 if (index0<0) index0=0;
2328 if (index1>=npoints-1) index1=npoints-1;
2329 if ((index1-index0)>1){
2330 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2331 }else{
2332 outy[ipoint]=graph->GetY()[ipoint];
2333 }
2334 }
0d1b4cf8 2335 // TLinearFitter fitter(3,"pol2");
2336// for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2337// Double_t lx=graph->GetX()[ipoint];
2338// Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2339// Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2340// if (index0<0) index0=0;
2341// if (index1>=npoints-1) index1=npoints-1;
2342// fitter.ClearPoints();
2343// for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2344// if ((index1-index0)>1){
2345// outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2346// }else{
2347// outy[ipoint]=graph->GetY()[ipoint];
2348// }
2349// }
2350
2351
2352
a23ba1c3 2353 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2354 graph->GetY()[ipoint] = outy[ipoint];
2355 }
2356 delete[] outy;
2357}
2358
8166d5d7 2359Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
a23ba1c3 2360 //
2361 // Use constant interpolation outside of range
2362 //
2363 if (!graph) {
2364 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2365 return 0;
2366 }
5647625c 2367
a23ba1c3 2368 if (graph->GetN()<1){
5647625c 2369 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
a23ba1c3 2370 return 0;
2371 }
5647625c 2372
2373
a23ba1c3 2374 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2375 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
5647625c 2376
2b68ab9c 2377 // printf("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref));
5647625c 2378
2379 if(graph->GetN()==1)
2380 return graph->Eval(graph->GetX()[0]);
2381
2382
2383 return graph->Eval(xref);
a23ba1c3 2384}
2385
108a67db 2386Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2387 //
2388 // Use constant interpolation outside of range also for spline fits
2389 //
2390 if (!graph) {
2391 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2392 return 0;
2393 }
2394 if (graph->GetKnots()<1){
2395 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2396 return 0;
2397 }
2398 if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2399 if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2400 return graph->Eval( xref);
2401}
2402
1e722a63 2403Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2404 //
2405 // Filter DCS sensor information
2406 // ymin - minimal value
2407 // ymax - max value
2408 // maxdy - maximal deirivative
2409 // sigmaCut - cut on values and derivative in terms of RMS distribution
2410 // Return value - accepted fraction
2411 //
2412 // Algorithm:
2413 //
2414 // 0. Calculate median and rms of values in specified range
2415 // 1. Filter out outliers - median+-sigmaCut*rms
2416 // values replaced by median
2417 //
2418 AliSplineFit * fit = sensor->GetFit();
2419 if (!fit) return 0.;
2420 Int_t nknots = fit->GetKnots();
2421 if (nknots==0) {
2422 delete fit;
2423 sensor->SetFit(0);
2424 return 0;
2425 }
2426 //
2427 Double_t *yin0 = new Double_t[nknots];
2428 Double_t *yin1 = new Double_t[nknots];
2429 Int_t naccept=0;
2430
2431 for (Int_t iknot=0; iknot< nknots; iknot++){
2432 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2433 yin0[naccept] = fit->GetY0()[iknot];
2434 yin1[naccept] = fit->GetY1()[iknot];
2435 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2436 naccept++;
2437 }
2438 }
2439 if (naccept<1) {
2440 delete fit;
2441 sensor->SetFit(0);
a9f487a1 2442 delete [] yin0;
2443 delete [] yin1;
1e722a63 2444 return 0.;
2445 }
0d1b4cf8 2446
1e722a63 2447 Double_t medianY0=0, medianY1=0;
2448 Double_t rmsY0 =0, rmsY1=0;
2449 medianY0 = TMath::Median(naccept, yin0);
2450 medianY1 = TMath::Median(naccept, yin1);
2451 rmsY0 = TMath::RMS(naccept, yin0);
2452 rmsY1 = TMath::RMS(naccept, yin1);
2453 naccept=0;
2454 //
2455 // 1. Filter out outliers - median+-sigmaCut*rms
2456 // values replaced by median
2457 // if replaced the derivative set to 0
2458 //
2459 for (Int_t iknot=0; iknot< nknots; iknot++){
2460 Bool_t isOK=kTRUE;
2461 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2462 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2463 if (nknots<2) fit->GetY1()[iknot]=0;
2464 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2465 if (!isOK){
2466 fit->GetY0()[iknot]=medianY0;
2467 fit->GetY1()[iknot]=0;
2468 }else{
2469 naccept++;
2470 }
2471 }
2472 delete [] yin0;
2473 delete [] yin1;
2474 return Float_t(naccept)/Float_t(nknots);
2475}
2476
2477Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2478 //
2479 // Filter temperature array
2480 // tempArray - array of temperatures -
2481 // ymin - minimal accepted temperature - default 15
2482 // ymax - maximal accepted temperature - default 22
2483 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2484 // return value - fraction of filtered sensors
2485 const Double_t kMaxDy=0.1;
2486 Int_t nsensors=tempArray->NumSensors();
2487 if (nsensors==0) return 0.;
2488 Int_t naccept=0;
2489 for (Int_t isensor=0; isensor<nsensors; isensor++){
2490 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2491 if (!sensor) continue;
2492 //printf("%d\n",isensor);
2493 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2494 if (sensor->GetFit()==0){
0d1b4cf8 2495 //delete sensor;
1e722a63 2496 tempArray->RemoveSensorNum(isensor);
2497 }else{
2498 naccept++;
2499 }
2500 }
2501 return Float_t(naccept)/Float_t(nsensors);
2502}
2503
a23ba1c3 2504
8166d5d7 2505void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2506 //
2507 // Filter CE data
1e722a63 2508 // Input parameters:
2509 // deltaT - smoothing window (in seconds)
2510 // cutAbs - max distance of the time info to the median (in time bins)
2511 // cutSigma - max distance (in the RMS)
2512 // pcstream - optional debug streamer to store original and filtered info
2513 // Hardwired parameters:
2514 // kMinPoints =10; // minimal number of points to define the CE
2515 // kMinSectors=12; // minimal number of sectors to define sideCE
2516 // Algorithm:
2517 // 0. Filter almost emty graphs (kMinPoints=10)
2518 // 1. calculate median and RMS per side
2519 // 2. Filter graphs - in respect with side medians
2520 // - cutAbs and cutDelta used
2521 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2522 // 4. Calculate mean for A side and C side
2523 //
2524 const Int_t kMinPoints =10; // minimal number of points to define the CE
2525 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2526 const Int_t kMinTime =400; // minimal arrival time of CE
2527 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
a23ba1c3 2528 Double_t medianY=0;
1e722a63 2529 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
a23ba1c3 2530 if (!cearray) return;
1e722a63 2531 Double_t tmin=-1;
2532 Double_t tmax=-1;
2533 //
2534 //
a23ba1c3 2535 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2cd3a697 2536 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
a23ba1c3 2537 if ( tempMapCE && cavernPressureCE){
1e722a63 2538 //
c5dfa2dc 2539 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2540 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2541 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2542 Bool_t isOK=kTRUE;
1e722a63 2543 if (isOK) {
2544 // recalculate P/T correction map for time of the CE
2545 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2546 driftCalib->SetName("driftPTCE");
2547 driftCalib->SetTitle("driftPTCE");
2548 cearray->AddLast(driftCalib);
2549 }
a23ba1c3 2550 }
1e722a63 2551 //
2552 // 0. Filter almost emty graphs
2553 //
a23ba1c3 2554
1e722a63 2555 for (Int_t i=0; i<72;i++){
a23ba1c3 2556 TGraph *graph= (TGraph*)arrT->At(i);
29797f44 2557 if (!graph) continue;
2558 graph->Sort();
a23ba1c3 2559 if (graph->GetN()<kMinPoints){
2560 arrT->AddAt(0,i);
2561 delete graph; // delete empty graph
2562 continue;
2563 }
1e722a63 2564 if (tmin<0) tmin = graph->GetX()[0];
2565 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2566 //
2567 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2568 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2569 }
2570 //
2571 // 1. calculate median and RMS per side
2572 //
2573 TArrayF arrA(100000), arrC(100000);
2574 Int_t nA=0, nC=0;
2575 Double_t medianA=0, medianC=0;
2576 Double_t rmsA=0, rmsC=0;
2577 for (Int_t isec=0; isec<72;isec++){
2578 TGraph *graph= (TGraph*)arrT->At(isec);
2579 if (!graph) continue;
2580 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2581 if (graph->GetY()[ipoint]<kMinTime) continue;
2582 if (nA>=arrA.fN) arrA.Set(nA*2);
2583 if (nC>=arrC.fN) arrC.Set(nC*2);
2584 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2585 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2586 }
2587 }
2588 if (nA>0){
2589 medianA=TMath::Median(nA,arrA.fArray);
2590 rmsA =TMath::RMS(nA,arrA.fArray);
2591 }
2592 if (nC>0){
2593 medianC=TMath::Median(nC,arrC.fArray);
2594 rmsC =TMath::RMS(nC,arrC.fArray);
2595 }
2596 //
2597 // 2. Filter graphs - in respect with side medians
2598 //
2599 TArrayD vecX(100000), vecY(100000);
2600 for (Int_t isec=0; isec<72;isec++){
2601 TGraph *graph= (TGraph*)arrT->At(isec);
2602 if (!graph) continue;
2603 Double_t median = (isec%36<18) ? medianA: medianC;
2604 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2605 Int_t naccept=0;
29797f44 2606 // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2607 for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
1e722a63 2608 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2609 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2610 vecX[naccept]= graph->GetX()[ipoint];
2611 vecY[naccept]= graph->GetY()[ipoint];
2612 naccept++;
2613 }
2614 if (naccept<kMinPoints){
2615 arrT->AddAt(0,isec);
2616 delete graph; // delete empty graph
2617 continue;
2618 }
2619 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2620 delete graph;
2621 arrT->AddAt(graph2,isec);
2622 }
2623 //
2624 // 3. Cut in respect wit the graph median
2625 //
2626 for (Int_t i=0; i<72;i++){
2627 TGraph *graph= (TGraph*)arrT->At(i);
2628 if (!graph) continue;
2629 //
2630 // filter in range
2631 //
2632 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
a23ba1c3 2633 if (!graphTS0) continue;
2634 if (graphTS0->GetN()<kMinPoints) {
2635 delete graphTS0;
2636 delete graph;
2637 arrT->AddAt(0,i);
2638 continue;
2639 }
1e722a63 2640 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
f03e487b 2641 if (!graphTS) continue;
a23ba1c3 2642 graphTS->Sort();
2643 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2644 if (pcstream){
2645 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2646 (*pcstream)<<"filterCE"<<
2647 "run="<<run<<
2648 "isec="<<i<<
2649 "mY="<<medianY<<
2650 "graph.="<<graph<<
2651 "graphTS0.="<<graphTS0<<
2652 "graphTS.="<<graphTS<<
2653 "\n";
2654 }
2655 delete graphTS0;
a23ba1c3 2656 arrT->AddAt(graphTS,i);
2657 delete graph;
2658 }
1e722a63 2659 //
2660 // Recalculate the mean time A side C side
2661 //
2662 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2663 Int_t meanPoints=(nA+nC)/72; // mean number of points
2664 for (Int_t itime=0; itime<200; itime++){
2665 nA=0, nC=0;
2666 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2667 for (Int_t i=0; i<72;i++){
2668 TGraph *graph= (TGraph*)arrT->At(i);
2669 if (!graph) continue;
2670 if (graph->GetN()<(meanPoints/4)) continue;
2671 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2672 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2673 }
2674 xA[itime]=time;
2675 xC[itime]=time;
2676 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2677 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2678 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2679 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2680 }
2681 //
2682 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2683 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2684 if (pcstream){
2685 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2686 (*pcstream)<<"filterAC"<<
2687 "run="<<run<<
2688 "nA="<<nA<<
2689 "nC="<<nC<<
2690 "rmsTA="<<rmsTA<<
2691 "rmsTC="<<rmsTC<<
2692 "\n";
2693 }
2694 //
2695 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2696 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2697 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2698 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2699 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2700 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2701 delete grA;
2702 delete grC;
2703 if (nA<kMinSectors) arrT->AddAt(0,72);
2704 else arrT->AddAt(graphTSA,72);
2705 if (nC<kMinSectors) arrT->AddAt(0,73);
2706 else arrT->AddAt(graphTSC,73);
a23ba1c3 2707}
2708
2709
8166d5d7 2710void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2711 //
2712 // Filter Drift velocity measurement using the tracks
2713 // 0. remove outlyers - error based
2714 // cutSigma
2715 //
2716 //
2717 const Int_t kMinPoints=1; // minimal number of points to define value
2718 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2719 Double_t medianY=0;
2720 if (!arrT) return;
2721 for (Int_t i=0; i<arrT->GetEntries();i++){
14301155 2722 TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
a23ba1c3 2723 if (!graph) continue;
2724 if (graph->GetN()<kMinPoints){
2725 delete graph;
2726 arrT->AddAt(0,i);
2727 continue;
2728 }
5647625c 2729 TGraphErrors *graph2 = NULL;
2730 if(graph->GetN()<10) {
2731 graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2732 if (!graph2) {
2733 delete graph; arrT->AddAt(0,i); continue;
2734 }
2735 }
2736 else {
2737 graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2738 if (!graph2) {
2739 delete graph; arrT->AddAt(0,i); continue;
2740 }
a23ba1c3 2741 }
2742 if (graph2->GetN()<1) {
2743 delete graph; arrT->AddAt(0,i); continue;
2744 }
2745 graph2->SetName(graph->GetName());
2746 graph2->SetTitle(graph->GetTitle());
2747 arrT->AddAt(graph2,i);
2748 if (pcstream){
2749 (*pcstream)<<"filterTracks"<<
2750 "run="<<run<<
2751 "isec="<<i<<
2752 "mY="<<medianY<<
2753 "graph.="<<graph<<
2754 "graph2.="<<graph2<<
2755 "\n";
2756 }
2757 delete graph;
2758 }
2759}
2760
2761
a23ba1c3 2762
2763
2764
2765Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2766 //
2767 //
2768 // get laser time offset
2769 // median around timeStamp+-deltaT
2770 // QA - chi2 needed for later usage - to be added
2771 // - currently cut on error
2772 //
2773 Int_t kMinPoints=1;
2774 Double_t kMinDelay=0.01;
2775 Double_t kMinDelayErr=0.0001;
2776
2777 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2778 if (!array) return 0;
2779 TGraphErrors *tlaser=0;
2780 if (array){
2781 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2782 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2783 }
2784 if (!tlaser) return 0;
2785 Int_t npoints0= tlaser->GetN();
2786 if (npoints0==0) return 0;
2787 Double_t *xlaser = new Double_t[npoints0];
2788 Double_t *ylaser = new Double_t[npoints0];
2789 Int_t npoints=0;
2790 for (Int_t i=0;i<npoints0;i++){
2791 //printf("%d\n",i);
2792 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2793 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2794 xlaser[npoints]=tlaser->GetX()[npoints];
2795 ylaser[npoints]=tlaser->GetY()[npoints];
2796 npoints++;
2797 }
2798 //
2799 //
2800 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2801 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2802 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2803 if (index0<0) index0=0;
2804 if (index1>=npoints-1) index1=npoints-1;
f03e487b 2805 if (index1-index0<kMinPoints) {
2806 delete [] ylaser;
2807 delete [] xlaser;
2808 return 0;
2809 }
a23ba1c3 2810 //
2811 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2812 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2813 delete [] ylaser;
2814 delete [] xlaser;
2815 return mean;
2816}
0d1b4cf8 2817
2818
2819
2820
8166d5d7 2821void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
0d1b4cf8 2822 //
2823 // Filter Goofie data
1fabc823 2824 // goofieArray - points will be filtered
2825 // deltaT - smmothing time window
2826 // cutSigma - outler sigma cut in rms
2827 // minVn, maxVd- range absolute cut for variable vd/pt
2828 // - to be tuned
0d1b4cf8 2829 //
0d1b4cf8 2830 // Ignore goofie if not enough points
2831 //
2832 const Int_t kMinPoints = 3;
2833 //
2834
2835 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2836 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2837 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2838 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2839 if (!graphvd) return;
2840 if (graphvd->GetN()<kMinPoints){
2841 delete graphvd;
2842 goofieArray->GetSensorNum(2)->SetGraph(0);
2843 return;
2844 }
2845 //
2846 // 1. Caluclate medians of critical variables
2847 // drift velcocity
2848 // P/T
2849 // area near peak
2850 // area far peak
2851 //
2852 Double_t medianpt=0;
1fabc823 2853 Double_t medianvd=0, sigmavd=0;
0d1b4cf8 2854 Double_t medianan=0;
2855 Double_t medianaf=0;
1fabc823 2856 Int_t entries=graphvd->GetN();
2857 Double_t yvdn[10000];
2858 Int_t nvd=0;
2859 //
2860 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2861 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2862 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2863 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2864 yvdn[nvd++]=graphvd->GetY()[ipoint];
2865 }
2866 if (nvd<kMinPoints){
2867 delete graphvd;
2868 goofieArray->GetSensorNum(2)->SetGraph(0);
2869 return;
2870 }
2871 //
2872 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2873 if (nuni>=kMinPoints){
2874 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2875 }else{
2876 medianvd = TMath::Median(nvd, yvdn);
2877 }
2878
0d1b4cf8 2879 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2880 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2881 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2882 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2883 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2884 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
0d1b4cf8 2885 delete graphpt0;
2886 delete graphpt1;
2887 delete graphan0;
2888 delete graphan1;
2889 delete graphaf0;
2890 delete graphaf1;
2891 //
2892 // 2. Make outlyer graph
2893 //
1fabc823 2894 Int_t nOK=0;
0d1b4cf8 2895 TGraph graphOut(*graphvd);
2896 for (Int_t i=0; i<entries;i++){
2897 //
2898 Bool_t isOut=kFALSE;
1fabc823 2899 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2900 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2901
2902 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2903 isOut|=kTRUE;
0d1b4cf8 2904 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
1fabc823 2905 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2906 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
0d1b4cf8 2907 graphOut.GetY()[i]= (isOut)?1:0;
1fabc823 2908 if (!isOut) nOK++;
0d1b4cf8 2909 }
1fabc823 2910 if (nOK<kMinPoints) {
2911 delete graphvd;
2912 goofieArray->GetSensorNum(2)->SetGraph(0);
2913 return;
2914 }
0d1b4cf8 2915 //
2916 // 3. Filter out outlyers - and smooth
2917 //
2918 TVectorF vmedianArray(goofieArray->NumSensors());
2919 TVectorF vrmsArray(goofieArray->NumSensors());
2920 Double_t xnew[10000];
2921 Double_t ynew[10000];
2922 TObjArray junk;
2923 junk.SetOwner(kTRUE);
2924 Bool_t isOK=kTRUE;
2925 //
2926 //
2927 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2928 isOK=kTRUE;
2929 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2930 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2931 //
2932 if (!sensor) continue;
2933 graphOld = sensor->GetGraph();
2934 if (graphOld) {
2935 sensor->SetGraph(0);
2936 Int_t nused=0;
2937 for (Int_t i=0;i<entries;i++){
2938 if (graphOut.GetY()[i]>0.5) continue;
2939 xnew[nused]=graphOld->GetX()[i];
2940 ynew[nused]=graphOld->GetY()[i];
2941 nused++;
2942 }
2943 graphNew = new TGraph(nused,xnew,ynew);
2944 junk.AddLast(graphNew);
2945 junk.AddLast(graphOld);
2946 Double_t median=0;
2947 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2948 if (graphNew0!=0){
2949 junk.AddLast(graphNew0);
2950 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2951 if (graphNew1!=0){
1fabc823 2952 junk.AddLast(graphNew1);
0d1b4cf8 2953 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2954 if (graphNew2!=0) {
1fabc823 2955 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
0d1b4cf8 2956 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2957 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2958 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
1fabc823 2959 printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
0d1b4cf8 2960 vmedianArray[isensor]=median;
0d1b4cf8 2961 //
2962 }
2963 }
2964 }
2965 }
2966 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2967 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2968 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2969 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2970 (*pcstream)<<"goofieA"<<
2971 Form("isOK_%d.=",isensor)<<isOK<<
2972 Form("s_%d.=",isensor)<<sensor<<
2973 Form("gr_%d.=",isensor)<<graphOld<<
2974 Form("gr0_%d.=",isensor)<<graphNew0<<
2975 Form("gr1_%d.=",isensor)<<graphNew1<<
2976 Form("gr2_%d.=",isensor)<<graphNew2;
1fabc823 2977 if (isOK) sensor->SetGraph(graphNew2);
2978 }
2979 (*pcstream)<<"goofieA"<<
2980 "vmed.="<<&vmedianArray<<
2981 "vrms.="<<&vrmsArray<<
2982 "\n";
2983 junk.Delete(); // delete temoprary graphs
0d1b4cf8 2984
2985}
2986
2987
a980538f 2988
2989
2990
8166d5d7 2991TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
a980538f 2992 //
2993 // Make a statistic matrix
2994 // Input parameters:
2995 // array - TObjArray of AliRelKalmanAlign
2996 // minFraction - minimal ration of accepted tracks
2997 // minStat - minimal statistic (number of accepted tracks)
2998 // maxvd - maximal deviation for the 1
2999 // Output matrix:
3000 // columns - Mean, Median, RMS
3001 // row - parameter type (rotation[3], translation[3], drift[3])
3002 if (!array) return 0;
3003 if (array->GetEntries()<=0) return 0;
cc65e4f5 3004 // Int_t entries = array->GetEntries();
a980538f 3005 Int_t entriesFast = array->GetEntriesFast();
3006 TVectorD state(9);
3007 TVectorD *valArray[9];
3008 for (Int_t i=0; i<9; i++){
3009 valArray[i] = new TVectorD(entriesFast);
3010 }
3011 Int_t naccept=0;
3012 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3013 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3014 if (!kalman) continue;
3015 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3016 if (kalman->GetNUpdates()<minStat) continue;
3017 if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
3018 kalman->GetState(state);
3019 for (Int_t ipar=0; ipar<9; ipar++)
3020 (*valArray[ipar])[naccept]=state[ipar];
3021 naccept++;
3022 }
5647625c 3023 //if (naccept<2) return 0;
3024 if (naccept<1) return 0;
a980538f 3025 TMatrixD *pstat=new TMatrixD(9,3);
3026 TMatrixD &stat=*pstat;
3027 for (Int_t ipar=0; ipar<9; ipar++){
3028 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3029 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3030 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3031 }
3032 return pstat;
3033}
3034
3035
8166d5d7 3036TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
a980538f 3037 //
3038 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3039 // Input:
3040 // array - input array
3041 // stat - mean parameters statistic
3042 // direction -
3043 // sigmaCut - maximal allowed deviation from mean in terms of RMS
3044 if (!array) return 0;
3045 if (array->GetEntries()<=0) return 0;
32100b1d 3046 if (!(&stat)) return 0;
cc65e4f5 3047 // error increase in 1 hour
3048 const Double_t kerrsTime[9]={
3049 0.00001, 0.00001, 0.00001,
3050 0.001, 0.001, 0.001,
dfb83639 3051 0.002, 0.01, 0.001};
cc65e4f5 3052 //
3053 //
a980538f 3054 Int_t entries = array->GetEntriesFast();
3055 TObjArray *sArray= new TObjArray(entries);
3056 AliRelAlignerKalman * sKalman =0;
3057 TVectorD state(9);
3058 for (Int_t i=0; i<entries; i++){
3059 Int_t index=(direction)? entries-i-1:i;
3060 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3061 if (!kalman) continue;
3062 Bool_t isOK=kTRUE;
3063 kalman->GetState(state);
3064 for (Int_t ipar=0; ipar<9; ipar++){
3065 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3066 }
3067 if (!sKalman &&isOK) {
3068 sKalman=new AliRelAlignerKalman(*kalman);
3069 sKalman->SetRejectOutliers(kFALSE);
3070 sKalman->SetRunNumber(kalman->GetRunNumber());
3071 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3072 }
3073 if (!sKalman) continue;
3074 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
cc65e4f5 3075 for (Int_t ipar=0; ipar<9; ipar++){
3076// (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3077// (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3078// (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3079 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3080 }
a980538f 3081 sKalman->SetRunNumber(kalman->GetRunNumber());
3082 if (!isOK) sKalman->SetRunNumber(0);
3083 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3084 if (!isOK) continue;
3085 sKalman->SetRejectOutliers(kFALSE);
3086 sKalman->SetRunNumber(kalman->GetRunNumber());
3087 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3088 sKalman->Merge(kalman);
3089 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3090 //sKalman->Print();
3091 }
3092 return sArray;
3093}
3094
8166d5d7 3095TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
cc65e4f5 3096 //
3097 // Merge 2 RelKalman arrays
3098 // Input:
3099 // arrayP - rel kalman in direction plus
3100 // arrayM - rel kalman in direction minus
3101 if (!arrayP) return 0;
3102 if (arrayP->GetEntries()<=0) return 0;
3103 if (!arrayM) return 0;
3104 if (arrayM->GetEntries()<=0) return 0;
3a4b0aed 3105
cc65e4f5 3106 Int_t entries = arrayP->GetEntriesFast();
3107 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3a4b0aed 3108
cc65e4f5 3109 for (Int_t i=0; i<entries; i++){
3110 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3111 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3112 if (!kalmanP) continue;
3113 if (!kalmanM) continue;
3a4b0aed 3114
3115 AliRelAlignerKalman *kalman = NULL;
3116 if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3117 kalman = new AliRelAlignerKalman(*kalmanP);
3118 kalman->Merge(kalmanM);
3119 }
3120 else if (kalmanP->GetRunNumber() == 0) {
3121 kalman = new AliRelAlignerKalman(*kalmanM);
3122 }
3123 else if (kalmanM->GetRunNumber() == 0) {
3124 kalman = new AliRelAlignerKalman(*kalmanP);
3125 }
3126 else
3127 continue;
3128
cc65e4f5 3129 array->AddAt(kalman,i);
3130 }
3a4b0aed 3131
cc65e4f5 3132 return array;
3133}