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