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