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