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