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