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