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