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