]> 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
Minor corrections to air-filled container volumes (M. Sitta)
[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 //
1411 if (!rocData) {
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);
1518 Int_t isOut=0;
1519 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1520 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1521 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1522 rocOut->SetValue(ichannel,isOut);
1523 if (isOut) noutliersPulser++;
1524 }
1525 }
1526 return out;
1527}
1528
1529
1530AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1531 //
1532 // Author : Marian Ivanov
1533 // Create pad time0 correction map using information from the CE and from pulser
1534 //
1535 //
1536 // Return PadTime0 to be used for time0 relative alignment
1537 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1538 // using $ALICE_ROOT/TPC/script/gui application
1539 //
1540 // fitResultsA - fitParameters A side
1541 // fitResultsC - fitParameters C side
1542 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1543 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1544 //
1545 //
1546 // Algorithm:
1547 // 1. Find outlier map for CE
1548 // 2. Find outlier map for Pulser
1549 // 3. Replace outlier by median at given sector (median without outliers)
1550 // 4. Substract from the CE data pulser
1551 // 5. Fit the CE with formula
1552 // 5.1) (IROC-OROC) offset
1553 // 5.2) gx
1554 // 5.3) gy
1555 // 5.4) (lx-xmid)
1556 // 5.5) (IROC-OROC)*(lx-xmid)
1557 // 5.6) (ly/lx)^2
1558 // 6. Substract gy fit dependence from the CE data
1559 // 7. Add pulser back to CE data
1560 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1561 // 9. return CE data
1562 //
1563 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1564 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1565
1566 // fit formula
1567 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)";
1568 // output for fit formula
1569 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)";
1570 // gy part of formula
1571 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)";
1572 //
1573 //
1574 if (!fCETmean) return 0;
1575 Double_t pgya,pgyc,pchi2a,pchi2c;
1576 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1577 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1578
1579 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1580 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1581 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1582 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1583 padPulser->SetName("padPulser");
1584 padPulserOut->SetName("padPulserOut");
1585 padCE->SetName("padCE");
1586 padCEIn->SetName("padCEIn");
1587 padCEOut->SetName("padCEOut");
1588 padOut->SetName("padOut");
1589
1590 //
1591 // make combined outlyers map
1592 // and replace outlyers in maps with median for chamber
1593 //
1594 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1595 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1596 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1597 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1598 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1599 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1600 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1601 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1602 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1603 if (rocPulserOut->GetValue(ichannel)>0) {
1604 rocPulser->SetValue(ichannel,pulserMedian);
1605 rocOut->SetValue(ichannel,1);
1606 }
1607 if (rocCEOut->GetValue(ichannel)>0) {
1608 rocCE->SetValue(ichannel,ceMedian);
1609 rocOut->SetValue(ichannel,1);
1610 }
1611 }
1612 }
1613 //
1614 // remove pulser time 0
1615 //
1616 padCE->Add(padPulser,-1);
1617 //
1618 // Make fits
1619 //
1620 TMatrixD dummy;
1621 Float_t chi2Af,chi2Cf;
1622 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1623 chi2A=chi2Af;
1624 chi2C=chi2Cf;
1625 //
1626 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1627 padCEFitGY->SetName("padCEFitGy");
1628 //
1629 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1630 padCEFit->SetName("padCEFit");
1631 //
1632 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1633 padCEDiff->SetName("padCEDiff");
1634 padCEDiff->Add(padCEFit,-1.);
1635 //
1636 //
1637 padCE->Add(padCEFitGY,-1.);
1638
1639 padCE->Add(padPulser,1.);
1640 Double_t padmedian = padCE->GetMedian();
1641 padCE->Add(-padmedian); // normalize to median
1642 //
1643 // Replace outliers by fit value - median of diff per given chamber -GY fit
1644 //
1645 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1646 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1647 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1648 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1649 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1650 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1651 //
1652 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1653 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1654 if (rocOut->GetValue(ichannel)==0) continue;
1655 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1656 rocCE->SetValue(ichannel,value);
1657 }
1658 }
1659 //
1660 //
1661 if (dumpfile){
1662 //dump to the file - result can be visualized
1663 AliTPCPreprocessorOnline preprocesor;
1664 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1665 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1666 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1667 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1668 //
1669 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1670 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1671 //
1672 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1673 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1674 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1675 preprocesor.DumpToFile(dumpfile);
1676 }
1677 delete padPulser;
1678 delete padPulserOut;
1679 delete padCEIn;
1680 delete padCEOut;
1681 delete padOut;
1682 delete padCEDiff;
1683 delete padCEFitGY;
1684 return padCE;
1685}
1686
817766d5 1687
1688
1689
1690
1691Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1692 //
1693 // find the closest point to xref in x direction
1694 // return dx and value
5647625c 1695 dx = 0;
1696 y = 0;
1697
1698 if(!graph) return 0;
1699 if(graph->GetN() < 1) return 0;
1700
817766d5 1701 Int_t index=0;
1702 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1703 if (index<0) index=0;
5647625c 1704 if(graph->GetN()==1) {
1705 dx = xref-graph->GetX()[index];
1706 }
1707 else {
1708 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1709 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1710 dx = xref-graph->GetX()[index];
1711 }
817766d5 1712 y = graph->GetY()[index];
1713 return index;
1714}
1715
817766d5 1716Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1717 //
1718 // Get the correction of the trigger offset
1719 // combining information from the laser track calibration
1720 // and from cosmic calibration
1721 //
1722 // run - run number
1723 // timeStamp - tim stamp in seconds
1724 // deltaT - integration period to calculate offset
1725 // deltaTLaser -max validity of laser data
1726 // valType - 0 - median, 1- mean
1727 //
1728 // Integration vaues are just recomendation - if not possible to get points
1729 // automatically increase the validity by factor 2
1730 // (recursive algorithm until one month of data taking)
1731 //
1732 //
1733 const Float_t kLaserCut=0.0005;
aec07590 1734 const Int_t kMaxPeriod=3600*24*30*12; // one year max
817766d5 1735 const Int_t kMinPoints=20;
1736 //
1737 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1738 if (!array) {
1739 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1740 }
1741 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1742 if (!array) return 0;
1743 //
1744 TGraphErrors *laserA[3]={0,0,0};
1745 TGraphErrors *laserC[3]={0,0,0};
1746 TGraphErrors *cosmicAll=0;
1747 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1748 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1749 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1750 //
1751 //
1752 if (!cosmicAll) return 0;
1753 Int_t nmeasC=cosmicAll->GetN();
1754 Float_t *tdelta = new Float_t[nmeasC];
1755 Int_t nused=0;
1756 for (Int_t i=0;i<nmeasC;i++){
1757 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1758 Float_t ccosmic=cosmicAll->GetY()[i];
1759 Double_t yA=0,yC=0,dA=0,dC=0;
1760 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1761 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1762 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1763 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1764 //
1765 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1766 Float_t claser=0;
1767 if (TMath::Abs(yA-yC)<kLaserCut) {
1768 claser=(yA-yC)*0.5;
1769 }else{
1770 if (i%2==0) claser=yA;
1771 if (i%2==1) claser=yC;
1772 }
1773 tdelta[nused]=ccosmic-claser;
1774 nused++;
1775 }
1776 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
aec07590 1777 if (nused<kMinPoints) {
1778 printf("AliFatal: No time offset calibration available\n");
1779 return 0;
1780 }
817766d5 1781 Double_t median = TMath::Median(nused,tdelta);
1782 Double_t mean = TMath::Mean(nused,tdelta);
a9f487a1 1783 delete [] tdelta;
817766d5 1784 return (valType==0) ? median:mean;
1785}
1786
a23ba1c3 1787Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
817766d5 1788 //
1789 // Get the correction of the drift velocity
1790 // combining information from the laser track calibration
1791 // and from cosmic calibration
1792 //
a23ba1c3 1793 // dist - return value - distance to closest point in graph
817766d5 1794 // run - run number
1795 // timeStamp - tim stamp in seconds
1796 // deltaT - integration period to calculate time0 offset
1797 // deltaTLaser -max validity of laser data
1798 // valType - 0 - median, 1- mean
1799 //
1800 // Integration vaues are just recomendation - if not possible to get points
1801 // automatically increase the validity by factor 2
1802 // (recursive algorithm until one month of data taking)
1803 //
1804 //
1805 //
1806 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1807 if (!array) {
1808 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1809 }
1810 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1811 if (!array) return 0;
1812 TGraphErrors *cosmicAll=0;
1813 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1814 if (!cosmicAll) return 0;
a23ba1c3 1815 Double_t grY=0;
1816 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1817
817766d5 1818 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
5647625c 1819 Double_t vcosmic = AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1e722a63 1820 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
817766d5 1821 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
cc65e4f5 1822 return vcosmic-t0;
817766d5 1823
1824 /*
1825 Example usage:
1826
1827 Int_t run=89000
1828 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1829 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1830 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1831 //
1832 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1833 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1834 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1835 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1836
1837 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1838 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1839
1840 */
1841
1842}
1843
949d8707 1844const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1845{
1846 //
1847 // Create a default name for the gui file
1848 //
1849
1850 return Form("guiRefTreeRun%s.root",GetRefValidity());
1851}
1852
1853Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1854{
1855 //
1856 // Create a gui reference tree
1857 // if dirname and filename are empty default values will be used
1858 // this is the recommended way of using this function
1859 // it allows to check whether a file with the given run validity alredy exists
1860 //
1861 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1862 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1863 return kFALSE;
1864 }
1865
1866 TString file=filename;
1867 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1868
1869 AliTPCPreprocessorOnline prep;
1870 //noise and pedestals
1871 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1872 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1873 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1874 //pulser data
1875 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1876 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1877 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1878 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1879 //CE data
1880 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1881 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1882 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1883 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1884 //Altro data
1885 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1886 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1887 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1888 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1889 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1890 //QA
1891 AliTPCdataQA *dataQA=fRefDataQA;
1892 if (dataQA) {
1893 if (dataQA->GetNLocalMaxima())
1894 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1895 if (dataQA->GetMaxCharge())
1896 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1897 if (dataQA->GetMeanCharge())
1898 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1899 if (dataQA->GetNoThreshold())
1900 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1901 if (dataQA->GetNTimeBins())
1902 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1903 if (dataQA->GetNPads())
1904 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1905 if (dataQA->GetTimePosition())
1906 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1907 }
1908 prep.DumpToFile(file.Data());
1909 return kTRUE;
1910}
1911
a23ba1c3 1912Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1913 //
78f17711 1914 // Get the correction of the drift velocity using the offline laser tracks calbration
a23ba1c3 1915 //
1916 // run - run number
1917 // timeStamp - tim stamp in seconds
1918 // deltaT - integration period to calculate time0 offset
1919 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1920 // Note in case no data form both A and C side - the value from active side used
1921 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
78f17711 1922
1923 return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1924}
1925
1926Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1927 //
1928 // Get the correction of the drift velocity using the online laser tracks calbration
1929 //
1930 // run - run number
1931 // timeStamp - tim stamp in seconds
1932 // deltaT - integration period to calculate time0 offset
1933 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1934 // Note in case no data form both A and C side - the value from active side used
1935 TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1936
1937 Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1938 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1939 if (!param) return 0;
1940
1941 //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1942 dv*=param->GetDriftV()/2.61301900000000000e+06;
1943 if (dv>1e-20) dv=1/dv-1;
1944 else return 0;
1945 // T/P correction
1946 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1947
1948 AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
2cd3a697 1949 AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
78f17711 1950
1951 Double_t corrPTA=0;
1952 Double_t corrPTC=0;
1953
1954 if (temp&&press) {
1955 AliTPCCalibVdrift corr(temp,press,0);
1956 corrPTA=corr.GetPTRelative(timeStamp,0);
1957 corrPTC=corr.GetPTRelative(timeStamp,1);
1958 }
1959
1960 if (side==0) dv -= corrPTA;
1961 if (side==1) dv -= corrPTC;
1962 if (side==2) dv -= (corrPTA+corrPTC)/2;
1963
1964 return dv;
1965}
1966
1967Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1968 Int_t side, TObjArray * const array){
1969 //
1970 // common drift velocity retrieval for online and offline method
1971 //
a23ba1c3 1972 TGraphErrors *grlaserA=0;
1973 TGraphErrors *grlaserC=0;
1974 Double_t vlaserA=0, vlaserC=0;
1975 if (!array) return 0;
1976 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1977 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1978 Double_t deltaY;
5647625c 1979 if (grlaserA && grlaserA->GetN()>0) {
a23ba1c3 1980 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1981 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1982 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1983 }
5647625c 1984 if (grlaserC && grlaserC->GetN()>0) {
a23ba1c3 1985 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1986 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1987 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1988 }
1989 if (side==0) return vlaserA;
1990 if (side==1) return vlaserC;
1991 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1992 if (!grlaserA) return vlaserC;
1993 if (!grlaserC) return vlaserA;
1994 return mdrift;
1995}
1996
1997
1998Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1999 //
2000 // Get the correction of the drift velocity using the CE laser data
2001 // combining information from the CE, laser track calibration
2002 // and P/T calibration
2003 //
2004 // run - run number
2005 // timeStamp - tim stamp in seconds
2006 // deltaT - integration period to calculate time0 offset
2007 // side - 0 - A side, 1 - C side, 2 - mean from both sides
2008 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
0d1b4cf8 2009 if (!arrT) return 0;
a23ba1c3 2010 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
2011 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2012 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2013 //
2014 //
2015 Double_t corrPTA = 0, corrPTC=0;
2016 Double_t ltime0A = 0, ltime0C=0;
2017 Double_t gry=0;
2018 Double_t corrA=0, corrC=0;
2019 Double_t timeA=0, timeC=0;
5647625c 2020 const Double_t kEpsilon = 0.00001;
a23ba1c3 2021 TGraph *graphA = (TGraph*)arrT->At(72);
2022 TGraph *graphC = (TGraph*)arrT->At(73);
2023 if (!graphA && !graphC) return 0.;
2024 if (graphA &&graphA->GetN()>0) {
2025 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2026 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2027 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2028 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2029 if(ltime0A < kEpsilon) return 0;
a23ba1c3 2030 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2031 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2032 corrA-=corrPTA;
2033 }
2034 if (graphC&&graphC->GetN()>0){
2035 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2036 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2037 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2038 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2039 if(ltime0C < kEpsilon) return 0;
2040if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2041 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2042 corrC-=corrPTC;
2043 }
2044
2045 if (side ==0 ) return corrA;
2046 if (side ==1 ) return corrC;
2047 Double_t corrM= (corrA+corrC)*0.5;
2048 if (!graphA) corrM=corrC;
2049 if (!graphC) corrM=corrA;
2050 return corrM;
2051}
2052
cc65e4f5 2053Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2054 //
2055 // return drift velocity using the TPC-ITS matchin method
2056 // return also distance to the closest point
2057 //
2058 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2059 TGraphErrors *graph=0;
2060 dist=0;
2061 if (!array) return 0;
5647625c 2062 //array->ls();
cc65e4f5 2063 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
2064 if (!graph) return 0;
2065 Double_t deltaY;
2066 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2067 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2068 return value;
2069}
2070
2071Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2072 //
2073 // Get time dependent time 0 (trigger delay in cm) correction
2074 // Arguments:
2075 // timestamp - timestamp
2076 // run - run number
2077 //
2078 // Notice - Extrapolation outside of calibration range - using constant function
2079 //
2080 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2081 TGraphErrors *graph=0;
2082 dist=0;
2083 if (!array) return 0;
2084 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
2085 if (!graph) return 0;
2086 Double_t deltaY;
2087 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2088 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2089 return value;
2090}
2091
2092
a23ba1c3 2093
2094
2095
2096Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2097 //
2098 // VERY obscure method - we need something in framework
2099 // Find the TPC runs with temperature OCDB entry
2100 // cache the start and end of the run
2101 //
2102 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2103 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2104 if (!storage) return 0;
2105 TString path=storage->GetURI();
2106 TString runsT;
2107 {
2108 TString command;
2109 if (path.Contains("local")){ // find the list if local system
2110 path.ReplaceAll("local://","");
2111 path+="TPC/Calib/Temperature";
2112 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2113 }
2114 runsT=gSystem->GetFromPipe(command);
2115 }
2116 TObjArray *arr= runsT.Tokenize("r");
2117 if (!arr) return 0;
2118 //
2119 TArrayI indexes(arr->GetEntries());
2120 TArrayI runs(arr->GetEntries());
2121 Int_t naccept=0;
2122 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2123 Int_t irunN = atoi(arr->At(irun)->GetName());
2124 if (irunN<startRun) continue;
2125 if (irunN>stopRun) continue;
2126 runs[naccept]=irunN;
2127 naccept++;
2128 }}
2129 fRuns.Set(naccept);
2130 fRunsStart.Set(fRuns.fN);
2131 fRunsStop.Set(fRuns.fN);
2132 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2133 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2134
2135 //
2136 AliCDBEntry * entry = 0;
2137 {for (Int_t irun=0;irun<fRuns.fN; irun++){
2138 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2139 if (!entry) continue;
2140 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2141 if (!tmpRun) continue;
2142 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2143 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2144 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
2145 }}
2146 return fRuns.fN;
2147}
2148
2149
2150Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2151 //
2152 // binary search - find the run for given time stamp
2153 //
2154 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2155 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2156 Int_t cindex = -1;
2157 for (Int_t index=index0; index<=index1; index++){
2158 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2159 if (debug) {
2160 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2161 }
2162 }
2163 if (cindex<0) cindex =(index0+index1)/2;
2164 if (cindex<0) {
2165 return 0;
2166 }
2167 return fRuns[cindex];
2168}
2169
2170
2171
2172
2173
2174TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2175 //
2176 // filter outlyer measurement
2177 // Only points around median +- sigmaCut filtered
2178 if (!graph) return 0;
2179 Int_t kMinPoints=2;
2180 Int_t npoints0 = graph->GetN();
2181 Int_t npoints=0;
2182 Float_t rmsY=0;
a23ba1c3 2183 //
2184 //
2185 if (npoints0<kMinPoints) return 0;
a9f487a1 2186
2187 Double_t *outx=new Double_t[npoints0];
2188 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2189 for (Int_t iter=0; iter<3; iter++){
2190 npoints=0;
2191 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2192 if (graph->GetY()[ipoint]==0) continue;
2193 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2194 outx[npoints] = graph->GetX()[ipoint];
2195 outy[npoints] = graph->GetY()[ipoint];
2196 npoints++;
2197 }
2198 if (npoints<=1) break;
2199 medianY =TMath::Median(npoints,outy);
2200 rmsY =TMath::RMS(npoints,outy);
2201 }
2202 TGraph *graphOut=0;
2203 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2204 delete [] outx;
2205 delete [] outy;
a23ba1c3 2206 return graphOut;
2207}
2208
2209
2210TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2211 //
2212 // filter outlyer measurement
2213 // Only points around median +- cut filtered
2214 if (!graph) return 0;
2215 Int_t kMinPoints=2;
2216 Int_t npoints0 = graph->GetN();
2217 Int_t npoints=0;
2218 Float_t rmsY=0;
a23ba1c3 2219 //
2220 //
2221 if (npoints0<kMinPoints) return 0;
a9f487a1 2222
2223 Double_t *outx=new Double_t[npoints0];
2224 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2225 for (Int_t iter=0; iter<3; iter++){
2226 npoints=0;
2227 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2228 if (graph->GetY()[ipoint]==0) continue;
2229 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2230 outx[npoints] = graph->GetX()[ipoint];
2231 outy[npoints] = graph->GetY()[ipoint];
2232 npoints++;
2233 }
2234 if (npoints<=1) break;
2235 medianY =TMath::Median(npoints,outy);
2236 rmsY =TMath::RMS(npoints,outy);
2237 }
2238 TGraph *graphOut=0;
2239 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2240 delete [] outx;
2241 delete [] outy;
a23ba1c3 2242 return graphOut;
2243}
2244
2245
2246
8166d5d7 2247TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
a23ba1c3 2248 //
2249 // filter outlyer measurement
2250 // Only points with normalized errors median +- sigmaCut filtered
2251 //
2252 Int_t kMinPoints=10;
2253 Int_t npoints0 = graph->GetN();
2254 Int_t npoints=0;
2255 Float_t medianErr=0, rmsErr=0;
a9f487a1 2256 //
2257 //
2258 if (npoints0<kMinPoints) return 0;
2259
a23ba1c3 2260 Double_t *outx=new Double_t[npoints0];
2261 Double_t *outy=new Double_t[npoints0];
2262 Double_t *erry=new Double_t[npoints0];
2263 Double_t *nerry=new Double_t[npoints0];
2264 Double_t *errx=new Double_t[npoints0];
a9f487a1 2265
a23ba1c3 2266 for (Int_t iter=0; iter<3; iter++){
2267 npoints=0;
2268 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2269 nerry[npoints] = graph->GetErrorY(ipoint);
2270 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2271 erry[npoints] = graph->GetErrorY(ipoint);
2272 outx[npoints] = graph->GetX()[ipoint];
2273 outy[npoints] = graph->GetY()[ipoint];
2274 errx[npoints] = graph->GetErrorY(ipoint);
2275 npoints++;
2276 }
2277 if (npoints==0) break;
2278 medianErr=TMath::Median(npoints,erry);
2279 medianY =TMath::Median(npoints,outy);
2280 rmsErr =TMath::RMS(npoints,erry);
2281 }
2282 TGraphErrors *graphOut=0;
2283 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2284 delete []outx;
2285 delete []outy;
a23ba1c3 2286 delete []erry;
a9f487a1 2287 delete []nerry;
2288 delete []errx;
a23ba1c3 2289 return graphOut;
2290}
2291
2292
2293void AliTPCcalibDButil::Sort(TGraph *graph){
2294 //
2295 // sort array - neccessay for approx
2296 //
2297 Int_t npoints = graph->GetN();
2298 Int_t *indexes=new Int_t[npoints];
2299 Double_t *outx=new Double_t[npoints];
2300 Double_t *outy=new Double_t[npoints];
2301 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2302 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2303 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2304 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2305 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2306
a9f487a1 2307 delete [] indexes;
2308 delete [] outx;
2309 delete [] outy;
a23ba1c3 2310}
2311void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2312 //
2313 // smmoth graph - mean on the interval
2314 //
2315 Sort(graph);
2316 Int_t npoints = graph->GetN();
2317 Double_t *outy=new Double_t[npoints];
0d1b4cf8 2318
a23ba1c3 2319 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2320 Double_t lx=graph->GetX()[ipoint];
2321 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2322 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2323 if (index0<0) index0=0;
2324 if (index1>=npoints-1) index1=npoints-1;
2325 if ((index1-index0)>1){
2326 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2327 }else{
2328 outy[ipoint]=graph->GetY()[ipoint];
2329 }
2330 }
0d1b4cf8 2331 // TLinearFitter fitter(3,"pol2");
2332// for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2333// Double_t lx=graph->GetX()[ipoint];
2334// Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2335// Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2336// if (index0<0) index0=0;
2337// if (index1>=npoints-1) index1=npoints-1;
2338// fitter.ClearPoints();
2339// for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2340// if ((index1-index0)>1){
2341// outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2342// }else{
2343// outy[ipoint]=graph->GetY()[ipoint];
2344// }
2345// }
2346
2347
2348
a23ba1c3 2349 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2350 graph->GetY()[ipoint] = outy[ipoint];
2351 }
2352 delete[] outy;
2353}
2354
8166d5d7 2355Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
a23ba1c3 2356 //
2357 // Use constant interpolation outside of range
2358 //
2359 if (!graph) {
2360 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2361 return 0;
2362 }
5647625c 2363
a23ba1c3 2364 if (graph->GetN()<1){
5647625c 2365 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
a23ba1c3 2366 return 0;
2367 }
5647625c 2368
2369
a23ba1c3 2370 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2371 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
5647625c 2372
2b68ab9c 2373 // printf("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref));
5647625c 2374
2375 if(graph->GetN()==1)
2376 return graph->Eval(graph->GetX()[0]);
2377
2378
2379 return graph->Eval(xref);
a23ba1c3 2380}
2381
108a67db 2382Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2383 //
2384 // Use constant interpolation outside of range also for spline fits
2385 //
2386 if (!graph) {
2387 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2388 return 0;
2389 }
2390 if (graph->GetKnots()<1){
2391 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2392 return 0;
2393 }
2394 if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2395 if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2396 return graph->Eval( xref);
2397}
2398
1e722a63 2399Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2400 //
2401 // Filter DCS sensor information
2402 // ymin - minimal value
2403 // ymax - max value
2404 // maxdy - maximal deirivative
2405 // sigmaCut - cut on values and derivative in terms of RMS distribution
2406 // Return value - accepted fraction
2407 //
2408 // Algorithm:
2409 //
2410 // 0. Calculate median and rms of values in specified range
2411 // 1. Filter out outliers - median+-sigmaCut*rms
2412 // values replaced by median
2413 //
2414 AliSplineFit * fit = sensor->GetFit();
2415 if (!fit) return 0.;
2416 Int_t nknots = fit->GetKnots();
2417 if (nknots==0) {
2418 delete fit;
2419 sensor->SetFit(0);
2420 return 0;
2421 }
2422 //
2423 Double_t *yin0 = new Double_t[nknots];
2424 Double_t *yin1 = new Double_t[nknots];
2425 Int_t naccept=0;
2426
2427 for (Int_t iknot=0; iknot< nknots; iknot++){
2428 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2429 yin0[naccept] = fit->GetY0()[iknot];
2430 yin1[naccept] = fit->GetY1()[iknot];
2431 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2432 naccept++;
2433 }
2434 }
2435 if (naccept<1) {
2436 delete fit;
2437 sensor->SetFit(0);
a9f487a1 2438 delete [] yin0;
2439 delete [] yin1;
1e722a63 2440 return 0.;
2441 }
0d1b4cf8 2442
1e722a63 2443 Double_t medianY0=0, medianY1=0;
2444 Double_t rmsY0 =0, rmsY1=0;
2445 medianY0 = TMath::Median(naccept, yin0);
2446 medianY1 = TMath::Median(naccept, yin1);
2447 rmsY0 = TMath::RMS(naccept, yin0);
2448 rmsY1 = TMath::RMS(naccept, yin1);
2449 naccept=0;
2450 //
2451 // 1. Filter out outliers - median+-sigmaCut*rms
2452 // values replaced by median
2453 // if replaced the derivative set to 0
2454 //
2455 for (Int_t iknot=0; iknot< nknots; iknot++){
2456 Bool_t isOK=kTRUE;
2457 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2458 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2459 if (nknots<2) fit->GetY1()[iknot]=0;
2460 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2461 if (!isOK){
2462 fit->GetY0()[iknot]=medianY0;
2463 fit->GetY1()[iknot]=0;
2464 }else{
2465 naccept++;
2466 }
2467 }
2468 delete [] yin0;
2469 delete [] yin1;
2470 return Float_t(naccept)/Float_t(nknots);
2471}
2472
2473Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2474 //
2475 // Filter temperature array
2476 // tempArray - array of temperatures -
2477 // ymin - minimal accepted temperature - default 15
2478 // ymax - maximal accepted temperature - default 22
2479 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2480 // return value - fraction of filtered sensors
2481 const Double_t kMaxDy=0.1;
2482 Int_t nsensors=tempArray->NumSensors();
2483 if (nsensors==0) return 0.;
2484 Int_t naccept=0;
2485 for (Int_t isensor=0; isensor<nsensors; isensor++){
2486 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2487 if (!sensor) continue;
2488 //printf("%d\n",isensor);
2489 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2490 if (sensor->GetFit()==0){
0d1b4cf8 2491 //delete sensor;
1e722a63 2492 tempArray->RemoveSensorNum(isensor);
2493 }else{
2494 naccept++;
2495 }
2496 }
2497 return Float_t(naccept)/Float_t(nsensors);
2498}
2499
a23ba1c3 2500
8166d5d7 2501void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2502 //
2503 // Filter CE data
1e722a63 2504 // Input parameters:
2505 // deltaT - smoothing window (in seconds)
2506 // cutAbs - max distance of the time info to the median (in time bins)
2507 // cutSigma - max distance (in the RMS)
2508 // pcstream - optional debug streamer to store original and filtered info
2509 // Hardwired parameters:
2510 // kMinPoints =10; // minimal number of points to define the CE
2511 // kMinSectors=12; // minimal number of sectors to define sideCE
2512 // Algorithm:
2513 // 0. Filter almost emty graphs (kMinPoints=10)
2514 // 1. calculate median and RMS per side
2515 // 2. Filter graphs - in respect with side medians
2516 // - cutAbs and cutDelta used
2517 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2518 // 4. Calculate mean for A side and C side
2519 //
2520 const Int_t kMinPoints =10; // minimal number of points to define the CE
2521 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2522 const Int_t kMinTime =400; // minimal arrival time of CE
2523 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
a23ba1c3 2524 Double_t medianY=0;
1e722a63 2525 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
a23ba1c3 2526 if (!cearray) return;
1e722a63 2527 Double_t tmin=-1;
2528 Double_t tmax=-1;
2529 //
2530 //
a23ba1c3 2531 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2cd3a697 2532 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
a23ba1c3 2533 if ( tempMapCE && cavernPressureCE){
1e722a63 2534 //
c5dfa2dc 2535 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2536 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2537 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2538 Bool_t isOK=kTRUE;
1e722a63 2539 if (isOK) {
2540 // recalculate P/T correction map for time of the CE
2541 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2542 driftCalib->SetName("driftPTCE");
2543 driftCalib->SetTitle("driftPTCE");
2544 cearray->AddLast(driftCalib);
2545 }
a23ba1c3 2546 }
1e722a63 2547 //
2548 // 0. Filter almost emty graphs
2549 //
a23ba1c3 2550
1e722a63 2551 for (Int_t i=0; i<72;i++){
a23ba1c3 2552 TGraph *graph= (TGraph*)arrT->At(i);
29797f44 2553 if (!graph) continue;
2554 graph->Sort();
a23ba1c3 2555 if (graph->GetN()<kMinPoints){
2556 arrT->AddAt(0,i);
2557 delete graph; // delete empty graph
2558 continue;
2559 }
1e722a63 2560 if (tmin<0) tmin = graph->GetX()[0];
2561 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2562 //
2563 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2564 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2565 }
2566 //
2567 // 1. calculate median and RMS per side
2568 //
2569 TArrayF arrA(100000), arrC(100000);
2570 Int_t nA=0, nC=0;
2571 Double_t medianA=0, medianC=0;
2572 Double_t rmsA=0, rmsC=0;
2573 for (Int_t isec=0; isec<72;isec++){
2574 TGraph *graph= (TGraph*)arrT->At(isec);
2575 if (!graph) continue;
2576 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2577 if (graph->GetY()[ipoint]<kMinTime) continue;
2578 if (nA>=arrA.fN) arrA.Set(nA*2);
2579 if (nC>=arrC.fN) arrC.Set(nC*2);
2580 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2581 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2582 }
2583 }
2584 if (nA>0){
2585 medianA=TMath::Median(nA,arrA.fArray);
2586 rmsA =TMath::RMS(nA,arrA.fArray);
2587 }
2588 if (nC>0){
2589 medianC=TMath::Median(nC,arrC.fArray);
2590 rmsC =TMath::RMS(nC,arrC.fArray);
2591 }
2592 //
2593 // 2. Filter graphs - in respect with side medians
2594 //
2595 TArrayD vecX(100000), vecY(100000);
2596 for (Int_t isec=0; isec<72;isec++){
2597 TGraph *graph= (TGraph*)arrT->At(isec);
2598 if (!graph) continue;
2599 Double_t median = (isec%36<18) ? medianA: medianC;
2600 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2601 Int_t naccept=0;
29797f44 2602 // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2603 for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
1e722a63 2604 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2605 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2606 vecX[naccept]= graph->GetX()[ipoint];
2607 vecY[naccept]= graph->GetY()[ipoint];
2608 naccept++;
2609 }
2610 if (naccept<kMinPoints){
2611 arrT->AddAt(0,isec);
2612 delete graph; // delete empty graph
2613 continue;
2614 }
2615 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2616 delete graph;
2617 arrT->AddAt(graph2,isec);
2618 }
2619 //
2620 // 3. Cut in respect wit the graph median
2621 //
2622 for (Int_t i=0; i<72;i++){
2623 TGraph *graph= (TGraph*)arrT->At(i);
2624 if (!graph) continue;
2625 //
2626 // filter in range
2627 //
2628 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
a23ba1c3 2629 if (!graphTS0) continue;
2630 if (graphTS0->GetN()<kMinPoints) {
2631 delete graphTS0;
2632 delete graph;
2633 arrT->AddAt(0,i);
2634 continue;
2635 }
1e722a63 2636 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
a23ba1c3 2637 graphTS->Sort();
2638 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2639 if (pcstream){
2640 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2641 (*pcstream)<<"filterCE"<<
2642 "run="<<run<<
2643 "isec="<<i<<
2644 "mY="<<medianY<<
2645 "graph.="<<graph<<
2646 "graphTS0.="<<graphTS0<<
2647 "graphTS.="<<graphTS<<
2648 "\n";
2649 }
2650 delete graphTS0;
2651 if (!graphTS) continue;
2652 arrT->AddAt(graphTS,i);
2653 delete graph;
2654 }
1e722a63 2655 //
2656 // Recalculate the mean time A side C side
2657 //
2658 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2659 Int_t meanPoints=(nA+nC)/72; // mean number of points
2660 for (Int_t itime=0; itime<200; itime++){
2661 nA=0, nC=0;
2662 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2663 for (Int_t i=0; i<72;i++){
2664 TGraph *graph= (TGraph*)arrT->At(i);
2665 if (!graph) continue;
2666 if (graph->GetN()<(meanPoints/4)) continue;
2667 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2668 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2669 }
2670 xA[itime]=time;
2671 xC[itime]=time;
2672 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2673 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2674 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2675 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2676 }
2677 //
2678 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2679 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2680 if (pcstream){
2681 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2682 (*pcstream)<<"filterAC"<<
2683 "run="<<run<<
2684 "nA="<<nA<<
2685 "nC="<<nC<<
2686 "rmsTA="<<rmsTA<<
2687 "rmsTC="<<rmsTC<<
2688 "\n";
2689 }
2690 //
2691 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2692 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2693 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2694 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2695 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2696 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2697 delete grA;
2698 delete grC;
2699 if (nA<kMinSectors) arrT->AddAt(0,72);
2700 else arrT->AddAt(graphTSA,72);
2701 if (nC<kMinSectors) arrT->AddAt(0,73);
2702 else arrT->AddAt(graphTSC,73);
a23ba1c3 2703}
2704
2705
8166d5d7 2706void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2707 //
2708 // Filter Drift velocity measurement using the tracks
2709 // 0. remove outlyers - error based
2710 // cutSigma
2711 //
2712 //
2713 const Int_t kMinPoints=1; // minimal number of points to define value
2714 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2715 Double_t medianY=0;
2716 if (!arrT) return;
2717 for (Int_t i=0; i<arrT->GetEntries();i++){
14301155 2718 TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
a23ba1c3 2719 if (!graph) continue;
2720 if (graph->GetN()<kMinPoints){
2721 delete graph;
2722 arrT->AddAt(0,i);
2723 continue;
2724 }
5647625c 2725 TGraphErrors *graph2 = NULL;
2726 if(graph->GetN()<10) {
2727 graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2728 if (!graph2) {
2729 delete graph; arrT->AddAt(0,i); continue;
2730 }
2731 }
2732 else {
2733 graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2734 if (!graph2) {
2735 delete graph; arrT->AddAt(0,i); continue;
2736 }
a23ba1c3 2737 }
2738 if (graph2->GetN()<1) {
2739 delete graph; arrT->AddAt(0,i); continue;
2740 }
2741 graph2->SetName(graph->GetName());
2742 graph2->SetTitle(graph->GetTitle());
2743 arrT->AddAt(graph2,i);
2744 if (pcstream){
2745 (*pcstream)<<"filterTracks"<<
2746 "run="<<run<<
2747 "isec="<<i<<
2748 "mY="<<medianY<<
2749 "graph.="<<graph<<
2750 "graph2.="<<graph2<<
2751 "\n";
2752 }
2753 delete graph;
2754 }
2755}
2756
2757
a23ba1c3 2758
2759
2760
2761Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2762 //
2763 //
2764 // get laser time offset
2765 // median around timeStamp+-deltaT
2766 // QA - chi2 needed for later usage - to be added
2767 // - currently cut on error
2768 //
2769 Int_t kMinPoints=1;
2770 Double_t kMinDelay=0.01;
2771 Double_t kMinDelayErr=0.0001;
2772
2773 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2774 if (!array) return 0;
2775 TGraphErrors *tlaser=0;
2776 if (array){
2777 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2778 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2779 }
2780 if (!tlaser) return 0;
2781 Int_t npoints0= tlaser->GetN();
2782 if (npoints0==0) return 0;
2783 Double_t *xlaser = new Double_t[npoints0];
2784 Double_t *ylaser = new Double_t[npoints0];
2785 Int_t npoints=0;
2786 for (Int_t i=0;i<npoints0;i++){
2787 //printf("%d\n",i);
2788 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2789 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2790 xlaser[npoints]=tlaser->GetX()[npoints];
2791 ylaser[npoints]=tlaser->GetY()[npoints];
2792 npoints++;
2793 }
2794 //
2795 //
2796 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2797 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2798 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2799 if (index0<0) index0=0;
2800 if (index1>=npoints-1) index1=npoints-1;
2801 if (index1-index0<kMinPoints) return 0;
2802 //
2803 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2804 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2805 delete [] ylaser;
2806 delete [] xlaser;
2807 return mean;
2808}
0d1b4cf8 2809
2810
2811
2812
8166d5d7 2813void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
0d1b4cf8 2814 //
2815 // Filter Goofie data
1fabc823 2816 // goofieArray - points will be filtered
2817 // deltaT - smmothing time window
2818 // cutSigma - outler sigma cut in rms
2819 // minVn, maxVd- range absolute cut for variable vd/pt
2820 // - to be tuned
0d1b4cf8 2821 //
0d1b4cf8 2822 // Ignore goofie if not enough points
2823 //
2824 const Int_t kMinPoints = 3;
2825 //
2826
2827 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2828 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2829 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2830 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2831 if (!graphvd) return;
2832 if (graphvd->GetN()<kMinPoints){
2833 delete graphvd;
2834 goofieArray->GetSensorNum(2)->SetGraph(0);
2835 return;
2836 }
2837 //
2838 // 1. Caluclate medians of critical variables
2839 // drift velcocity
2840 // P/T
2841 // area near peak
2842 // area far peak
2843 //
2844 Double_t medianpt=0;
1fabc823 2845 Double_t medianvd=0, sigmavd=0;
0d1b4cf8 2846 Double_t medianan=0;
2847 Double_t medianaf=0;
1fabc823 2848 Int_t entries=graphvd->GetN();
2849 Double_t yvdn[10000];
2850 Int_t nvd=0;
2851 //
2852 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2853 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2854 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2855 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2856 yvdn[nvd++]=graphvd->GetY()[ipoint];
2857 }
2858 if (nvd<kMinPoints){
2859 delete graphvd;
2860 goofieArray->GetSensorNum(2)->SetGraph(0);
2861 return;
2862 }
2863 //
2864 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2865 if (nuni>=kMinPoints){
2866 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2867 }else{
2868 medianvd = TMath::Median(nvd, yvdn);
2869 }
2870
0d1b4cf8 2871 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2872 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2873 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2874 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2875 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2876 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
0d1b4cf8 2877 delete graphpt0;
2878 delete graphpt1;
2879 delete graphan0;
2880 delete graphan1;
2881 delete graphaf0;
2882 delete graphaf1;
2883 //
2884 // 2. Make outlyer graph
2885 //
1fabc823 2886 Int_t nOK=0;
0d1b4cf8 2887 TGraph graphOut(*graphvd);
2888 for (Int_t i=0; i<entries;i++){
2889 //
2890 Bool_t isOut=kFALSE;
1fabc823 2891 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2892 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2893
2894 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2895 isOut|=kTRUE;
0d1b4cf8 2896 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
1fabc823 2897 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2898 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
0d1b4cf8 2899 graphOut.GetY()[i]= (isOut)?1:0;
1fabc823 2900 if (!isOut) nOK++;
0d1b4cf8 2901 }
1fabc823 2902 if (nOK<kMinPoints) {
2903 delete graphvd;
2904 goofieArray->GetSensorNum(2)->SetGraph(0);
2905 return;
2906 }
0d1b4cf8 2907 //
2908 // 3. Filter out outlyers - and smooth
2909 //
2910 TVectorF vmedianArray(goofieArray->NumSensors());
2911 TVectorF vrmsArray(goofieArray->NumSensors());
2912 Double_t xnew[10000];
2913 Double_t ynew[10000];
2914 TObjArray junk;
2915 junk.SetOwner(kTRUE);
2916 Bool_t isOK=kTRUE;
2917 //
2918 //
2919 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2920 isOK=kTRUE;
2921 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2922 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2923 //
2924 if (!sensor) continue;
2925 graphOld = sensor->GetGraph();
2926 if (graphOld) {
2927 sensor->SetGraph(0);
2928 Int_t nused=0;
2929 for (Int_t i=0;i<entries;i++){
2930 if (graphOut.GetY()[i]>0.5) continue;
2931 xnew[nused]=graphOld->GetX()[i];
2932 ynew[nused]=graphOld->GetY()[i];
2933 nused++;
2934 }
2935 graphNew = new TGraph(nused,xnew,ynew);
2936 junk.AddLast(graphNew);
2937 junk.AddLast(graphOld);
2938 Double_t median=0;
2939 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2940 if (graphNew0!=0){
2941 junk.AddLast(graphNew0);
2942 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2943 if (graphNew1!=0){
1fabc823 2944 junk.AddLast(graphNew1);
0d1b4cf8 2945 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2946 if (graphNew2!=0) {
1fabc823 2947 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
0d1b4cf8 2948 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2949 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2950 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
1fabc823 2951 printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
0d1b4cf8 2952 vmedianArray[isensor]=median;
0d1b4cf8 2953 //
2954 }
2955 }
2956 }
2957 }
2958 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2959 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2960 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2961 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2962 (*pcstream)<<"goofieA"<<
2963 Form("isOK_%d.=",isensor)<<isOK<<
2964 Form("s_%d.=",isensor)<<sensor<<
2965 Form("gr_%d.=",isensor)<<graphOld<<
2966 Form("gr0_%d.=",isensor)<<graphNew0<<
2967 Form("gr1_%d.=",isensor)<<graphNew1<<
2968 Form("gr2_%d.=",isensor)<<graphNew2;
1fabc823 2969 if (isOK) sensor->SetGraph(graphNew2);
2970 }
2971 (*pcstream)<<"goofieA"<<
2972 "vmed.="<<&vmedianArray<<
2973 "vrms.="<<&vrmsArray<<
2974 "\n";
2975 junk.Delete(); // delete temoprary graphs
0d1b4cf8 2976
2977}
2978
2979
a980538f 2980
2981
2982
8166d5d7 2983TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
a980538f 2984 //
2985 // Make a statistic matrix
2986 // Input parameters:
2987 // array - TObjArray of AliRelKalmanAlign
2988 // minFraction - minimal ration of accepted tracks
2989 // minStat - minimal statistic (number of accepted tracks)
2990 // maxvd - maximal deviation for the 1
2991 // Output matrix:
2992 // columns - Mean, Median, RMS
2993 // row - parameter type (rotation[3], translation[3], drift[3])
2994 if (!array) return 0;
2995 if (array->GetEntries()<=0) return 0;
cc65e4f5 2996 // Int_t entries = array->GetEntries();
a980538f 2997 Int_t entriesFast = array->GetEntriesFast();
2998 TVectorD state(9);
2999 TVectorD *valArray[9];
3000 for (Int_t i=0; i<9; i++){
3001 valArray[i] = new TVectorD(entriesFast);
3002 }
3003 Int_t naccept=0;
3004 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3005 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3006 if (!kalman) continue;
3007 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3008 if (kalman->GetNUpdates()<minStat) continue;
3009 if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
3010 kalman->GetState(state);
3011 for (Int_t ipar=0; ipar<9; ipar++)
3012 (*valArray[ipar])[naccept]=state[ipar];
3013 naccept++;
3014 }
5647625c 3015 //if (naccept<2) return 0;
3016 if (naccept<1) return 0;
a980538f 3017 TMatrixD *pstat=new TMatrixD(9,3);
3018 TMatrixD &stat=*pstat;
3019 for (Int_t ipar=0; ipar<9; ipar++){
3020 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3021 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3022 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3023 }
3024 return pstat;
3025}
3026
3027
8166d5d7 3028TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
a980538f 3029 //
3030 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3031 // Input:
3032 // array - input array
3033 // stat - mean parameters statistic
3034 // direction -
3035 // sigmaCut - maximal allowed deviation from mean in terms of RMS
3036 if (!array) return 0;
3037 if (array->GetEntries()<=0) return 0;
32100b1d 3038 if (!(&stat)) return 0;
cc65e4f5 3039 // error increase in 1 hour
3040 const Double_t kerrsTime[9]={
3041 0.00001, 0.00001, 0.00001,
3042 0.001, 0.001, 0.001,
dfb83639 3043 0.002, 0.01, 0.001};
cc65e4f5 3044 //
3045 //
a980538f 3046 Int_t entries = array->GetEntriesFast();
3047 TObjArray *sArray= new TObjArray(entries);
3048 AliRelAlignerKalman * sKalman =0;
3049 TVectorD state(9);
3050 for (Int_t i=0; i<entries; i++){
3051 Int_t index=(direction)? entries-i-1:i;
3052 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3053 if (!kalman) continue;
3054 Bool_t isOK=kTRUE;
3055 kalman->GetState(state);
3056 for (Int_t ipar=0; ipar<9; ipar++){
3057 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3058 }
3059 if (!sKalman &&isOK) {
3060 sKalman=new AliRelAlignerKalman(*kalman);
3061 sKalman->SetRejectOutliers(kFALSE);
3062 sKalman->SetRunNumber(kalman->GetRunNumber());
3063 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3064 }
3065 if (!sKalman) continue;
3066 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
cc65e4f5 3067 for (Int_t ipar=0; ipar<9; ipar++){
3068// (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3069// (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3070// (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3071 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3072 }
a980538f 3073 sKalman->SetRunNumber(kalman->GetRunNumber());
3074 if (!isOK) sKalman->SetRunNumber(0);
3075 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3076 if (!isOK) continue;
3077 sKalman->SetRejectOutliers(kFALSE);
3078 sKalman->SetRunNumber(kalman->GetRunNumber());
3079 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3080 sKalman->Merge(kalman);
3081 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3082 //sKalman->Print();
3083 }
3084 return sArray;
3085}
3086
8166d5d7 3087TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
cc65e4f5 3088 //
3089 // Merge 2 RelKalman arrays
3090 // Input:
3091 // arrayP - rel kalman in direction plus
3092 // arrayM - rel kalman in direction minus
3093 if (!arrayP) return 0;
3094 if (arrayP->GetEntries()<=0) return 0;
3095 if (!arrayM) return 0;
3096 if (arrayM->GetEntries()<=0) return 0;
3a4b0aed 3097
cc65e4f5 3098 Int_t entries = arrayP->GetEntriesFast();
3099 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3a4b0aed 3100
cc65e4f5 3101 for (Int_t i=0; i<entries; i++){
3102 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3103 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3104 if (!kalmanP) continue;
3105 if (!kalmanM) continue;
3a4b0aed 3106
3107 AliRelAlignerKalman *kalman = NULL;
3108 if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3109 kalman = new AliRelAlignerKalman(*kalmanP);
3110 kalman->Merge(kalmanM);
3111 }
3112 else if (kalmanP->GetRunNumber() == 0) {
3113 kalman = new AliRelAlignerKalman(*kalmanM);
3114 }
3115 else if (kalmanM->GetRunNumber() == 0) {
3116 kalman = new AliRelAlignerKalman(*kalmanP);
3117 }
3118 else
3119 continue;
3120
cc65e4f5 3121 array->AddAt(kalman,i);
3122 }
3a4b0aed 3123
cc65e4f5 3124 return array;
3125}