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