]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibDButil.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliTPCcalibDB - Applying filters for calration graphs using AliTPCcalibDButil
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibDButil.cxx
CommitLineData
892226be 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////////
18// //
19// Class providing the calculation of derived quantities (mean,rms,fits,...) //
20// of calibration entries //
21/*
22
23
24*/
25////////////////////////////////////////////////////////////////////////////////
26
27#include <TMath.h>
28#include <TVectorT.h>
29#include <TObjArray.h>
30#include <TGraph.h>
7390f655 31#include <TFile.h>
32#include <TDirectory.h>
a23ba1c3 33#include <AliCDBStorage.h>
892226be 34#include <AliDCSSensorArray.h>
a23ba1c3 35#include <AliTPCSensorTempArray.h>
892226be 36#include <AliDCSSensor.h>
37#include "AliTPCcalibDB.h"
38#include "AliTPCCalPad.h"
39#include "AliTPCCalROC.h"
40#include "AliTPCROC.h"
41#include "AliTPCmapper.h"
42#include "AliTPCParam.h"
6e7d7dc4 43#include "AliTPCCalibRaw.h"
817766d5 44#include "TGraphErrors.h"
a23ba1c3 45#include "AliLog.h"
892226be 46#include "AliTPCcalibDButil.h"
2cb269df 47#include "AliTPCPreprocessorOnline.h"
a23ba1c3 48#include "AliTPCCalibVdrift.h"
892226be 49
50ClassImp(AliTPCcalibDButil)
51AliTPCcalibDButil::AliTPCcalibDButil() :
52 TObject(),
1e722a63 53 fCalibDB(0),
892226be 54 fPadNoise(0x0),
55 fPedestals(0x0),
56 fPulserTmean(0x0),
57 fPulserTrms(0x0),
58 fPulserQmean(0x0),
59 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
60 fCETmean(0x0),
61 fCETrms(0x0),
62 fCEQmean(0x0),
63 fALTROMasked(0x0),
6e7d7dc4 64 fCalibRaw(0x0),
7390f655 65 fRefPadNoise(0x0),
66 fRefPedestals(0x0),
67 fRefPulserTmean(0x0),
68 fRefPulserTrms(0x0),
69 fRefPulserQmean(0x0),
70 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
71 fRefCETmean(0x0),
72 fRefCETrms(0x0),
73 fRefCEQmean(0x0),
74 fRefALTROMasked(0x0),
75 fRefCalibRaw(0x0),
892226be 76 fGoofieArray(0x0),
77 fMapper(new AliTPCmapper(0x0)),
78 fNpulserOutliers(-1),
7390f655 79 fIrocTimeOffset(0),
892226be 80 fCETmaxLimitAbs(1.5),
81 fPulTmaxLimitAbs(1.5),
82 fPulQmaxLimitAbs(5),
a23ba1c3 83 fPulQminLimit(11),
84 fRuns(0), // run list with OCDB info
85 fRunsStart(0), // start time for given run
86 fRunsStop(0) // stop time for given run
892226be 87{
88 //
89 // Default ctor
90 //
91}
92//_____________________________________________________________________________________
93AliTPCcalibDButil::~AliTPCcalibDButil()
94{
95 //
96 // dtor
97 //
98 delete fPulserOutlier;
7390f655 99 delete fRefPulserOutlier;
892226be 100 delete fMapper;
7390f655 101 if (fRefPadNoise) delete fRefPadNoise;
102 if (fRefPedestals) delete fRefPedestals;
103 if (fRefPulserTmean) delete fRefPulserTmean;
104 if (fRefPulserTrms) delete fRefPulserTrms;
105 if (fRefPulserQmean) delete fRefPulserQmean;
106 if (fRefCETmean) delete fRefCETmean;
107 if (fRefCETrms) delete fRefCETrms;
108 if (fRefCEQmean) delete fRefCEQmean;
109 if (fRefALTROMasked) delete fRefALTROMasked;
110 if (fRefCalibRaw) delete fRefCalibRaw;
111
892226be 112}
113//_____________________________________________________________________________________
114void AliTPCcalibDButil::UpdateFromCalibDB()
115{
116 //
117 // Update pointers from calibDB
118 //
1e722a63 119 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
892226be 120 fPadNoise=fCalibDB->GetPadNoise();
121 fPedestals=fCalibDB->GetPedestals();
122 fPulserTmean=fCalibDB->GetPulserTmean();
123 fPulserTrms=fCalibDB->GetPulserTrms();
124 fPulserQmean=fCalibDB->GetPulserQmean();
125 fCETmean=fCalibDB->GetCETmean();
126 fCETrms=fCalibDB->GetCETrms();
127 fCEQmean=fCalibDB->GetCEQmean();
128 fALTROMasked=fCalibDB->GetALTROMasked();
129 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
6e7d7dc4 130 fCalibRaw=fCalibDB->GetCalibRaw();
892226be 131 UpdatePulserOutlierMap();
132}
133//_____________________________________________________________________________________
7390f655 134void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
2cb269df 135 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad *outCE)
892226be 136{
137 //
138 // Process the CE data for this run
139 // the return TVectorD arrays contian the results of the fit
140 // noutliersCE contains the number of pads marked as outliers,
141 // not including masked and edge pads
142 //
143
144 //retrieve CE and ALTRO data
145 if (!fCETmean){
146 TString fitString(fitFormula);
147 fitString.ReplaceAll("++","#");
148 Int_t ndim=fitString.CountChar('#')+2;
149 fitResultsA.ResizeTo(ndim);
150 fitResultsC.ResizeTo(ndim);
151 fitResultsA.Zero();
152 fitResultsC.Zero();
153 noutliersCE=-1;
154 return;
155 }
156 noutliersCE=0;
157 //create outlier map
7390f655 158 AliTPCCalPad *out=0;
159 if (outCE) out=outCE;
160 else out=new AliTPCCalPad("outCE","outCE");
892226be 161 AliTPCCalROC *rocMasked=0x0;
162 //loop over all channels
163 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
164 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
165 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
7390f655 166 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
6e7d7dc4 167 if (!rocData) {
168 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
7390f655 169 rocOut->Add(1.);
6e7d7dc4 170 continue;
171 }
892226be 172 //add time offset to IROCs
173 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
174 rocData->Add(fIrocTimeOffset);
175 //select outliers
176 UInt_t nrows=rocData->GetNrows();
177 for (UInt_t irow=0;irow<nrows;++irow){
178 UInt_t npads=rocData->GetNPads(irow);
179 for (UInt_t ipad=0;ipad<npads;++ipad){
7390f655 180 rocOut->SetValue(irow,ipad,0);
892226be 181 //exclude masked pads
182 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
7390f655 183 rocOut->SetValue(irow,ipad,1);
892226be 184 continue;
185 }
732e90a8 186 //exclude first two rows in IROC and last two rows in OROC
187 if (iroc<36){
188 if (irow<2) rocOut->SetValue(irow,ipad,1);
189 } else {
190 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
191 }
892226be 192 //exclude edge pads
7390f655 193 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
194 Float_t valTmean=rocData->GetValue(irow,ipad);
892226be 195 //exclude values that are exactly 0
196 if (valTmean==0) {
7390f655 197 rocOut->SetValue(irow,ipad,1);
892226be 198 ++noutliersCE;
199 }
200 // exclude channels with too large variations
201 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
7390f655 202 rocOut->SetValue(irow,ipad,1);
892226be 203 ++noutliersCE;
204 }
205 }
206 }
207 }
208 //perform fit
209 TMatrixD dummy;
2cb269df 210 Float_t chi2Af,chi2Cf;
211 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
212 chi2A=chi2Af;
213 chi2C=chi2Cf;
7390f655 214 if (!outCE) delete out;
892226be 215}
216//_____________________________________________________________________________________
217void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
218 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
219 Float_t &driftTimeA, Float_t &driftTimeC )
220{
221 //
222 // Calculate statistical information from the CE graphs for drift time and charge
223 //
224
225 //reset arrays
226 vecTEntries.ResizeTo(72);
227 vecTMean.ResizeTo(72);
228 vecTRMS.ResizeTo(72);
229 vecTMedian.ResizeTo(72);
230 vecQEntries.ResizeTo(72);
231 vecQMean.ResizeTo(72);
232 vecQRMS.ResizeTo(72);
233 vecQMedian.ResizeTo(72);
234 vecTEntries.Zero();
235 vecTMean.Zero();
236 vecTRMS.Zero();
237 vecTMedian.Zero();
238 vecQEntries.Zero();
239 vecQMean.Zero();
240 vecQRMS.Zero();
241 vecQMedian.Zero();
242 driftTimeA=0;
243 driftTimeC=0;
244 TObjArray *arrT=fCalibDB->GetCErocTtime();
245 TObjArray *arrQ=fCalibDB->GetCErocQtime();
246 if (arrT){
247 for (Int_t isec=0;isec<74;++isec){
248 TGraph *gr=(TGraph*)arrT->At(isec);
249 if (!gr) continue;
250 TVectorD values;
251 Int_t npoints = gr->GetN();
252 values.ResizeTo(npoints);
253 Int_t nused =0;
6e7d7dc4 254 //skip first points, theres always some problems with finding the CE position
255 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
256 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
892226be 257 values[nused]=gr->GetY()[ipoint];
258 nused++;
259 }
260 }
261 //
262 if (isec<72) vecTEntries[isec]= nused;
263 if (nused>1){
264 if (isec<72){
265 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
266 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
267 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
268 } else if (isec==72){
269 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
270 } else if (isec==73){
271 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
272 }
273 }
274 }
275 }
276 if (arrQ){
6e7d7dc4 277 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
892226be 278 TGraph *gr=(TGraph*)arrQ->At(isec);
279 if (!gr) continue;
280 TVectorD values;
281 Int_t npoints = gr->GetN();
282 values.ResizeTo(npoints);
283 Int_t nused =0;
284 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
817766d5 285 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
892226be 286 values[nused]=gr->GetY()[ipoint];
287 nused++;
288 }
289 }
290 //
6e7d7dc4 291 vecQEntries[isec]= nused;
892226be 292 if (nused>1){
293 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
294 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
295 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
296 }
297 }
298 }
299}
300
301//_____________________________________________________________________________________
302void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
303 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
304 Int_t &nonMaskedZero)
305{
306 //
307 // process noise data
308 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
309 // OROCs small pads [2] and OROCs large pads [3]
310 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
311 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
312 //
313
314 //set proper size and reset
315 const UInt_t infoSize=4;
316 vNoiseMean.ResizeTo(infoSize);
317 vNoiseMeanSenRegions.ResizeTo(infoSize);
318 vNoiseRMS.ResizeTo(infoSize);
319 vNoiseRMSSenRegions.ResizeTo(infoSize);
320 vNoiseMean.Zero();
321 vNoiseMeanSenRegions.Zero();
322 vNoiseRMS.Zero();
323 vNoiseRMSSenRegions.Zero();
324 nonMaskedZero=0;
325 //counters
326 TVectorD c(infoSize);
327 TVectorD cs(infoSize);
328 //tpc parameters
329 AliTPCParam par;
330 par.Update();
331 //retrieve noise and ALTRO data
332 if (!fPadNoise) return;
333 AliTPCCalROC *rocMasked=0x0;
334 //create IROC, OROC1, OROC2 and sensitive region masks
335 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
336 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
337 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
338 UInt_t nrows=noiseROC->GetNrows();
339 for (UInt_t irow=0;irow<nrows;++irow){
340 UInt_t npads=noiseROC->GetNPads(irow);
341 for (UInt_t ipad=0;ipad<npads;++ipad){
342 //don't use masked channels;
343 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
344 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
345 //check if noise==0
346 if (noiseVal==0) {
347 ++nonMaskedZero;
348 continue;
349 }
350 //check for nan
351 if ( !(noiseVal<10000000) ){
352 printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
353 continue;
354 }
355 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
356 Int_t masksen=1; // sensitive pards are not masked (0)
357 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
358 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
359 //IROCs
360 if (irow>19&&irow<46){
361 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
362 }
363 Int_t type=1;
364 vNoiseMean[type]+=noiseVal;
365 vNoiseRMS[type]+=noiseVal*noiseVal;
366 ++c[type];
367 if (!masksen){
368 vNoiseMeanSenRegions[type]+=noiseVal;
369 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
370 ++cs[type];
371 }
372 } else {
373 //OROCs
374 //define sensive regions
375 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
376 if ( irow>75 ){
377 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
378 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
379 }
380 if ((Int_t)irow<par.GetNRowUp1()){
381 //OROC1
382 Int_t type=2;
383 vNoiseMean[type]+=noiseVal;
384 vNoiseRMS[type]+=noiseVal*noiseVal;
385 ++c[type];
386 if (!masksen){
387 vNoiseMeanSenRegions[type]+=noiseVal;
388 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
389 ++cs[type];
390 }
391 }else{
392 //OROC2
393 Int_t type=3;
394 vNoiseMean[type]+=noiseVal;
395 vNoiseRMS[type]+=noiseVal*noiseVal;
396 ++c[type];
397 if (!masksen){
398 vNoiseMeanSenRegions[type]+=noiseVal;
399 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
400 ++cs[type];
401 }
402 }
403 }
404 //whole tpc
405 Int_t type=0;
406 vNoiseMean[type]+=noiseVal;
407 vNoiseRMS[type]+=noiseVal*noiseVal;
408 ++c[type];
409 if (!masksen){
410 vNoiseMeanSenRegions[type]+=noiseVal;
411 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
412 ++cs[type];
413 }
414 }//end loop pads
415 }//end loop rows
416 }//end loop sectors (rocs)
417
418 //calculate mean and RMS
419 const Double_t verySmall=0.0000000001;
420 for (UInt_t i=0;i<infoSize;++i){
421 Double_t mean=0;
422 Double_t rms=0;
423 Double_t meanSen=0;
424 Double_t rmsSen=0;
425
426 if (c[i]>verySmall){
427// printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
428 mean=vNoiseMean[i]/c[i];
429 rms=vNoiseRMS[i];
430 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
431 }
432 vNoiseMean[i]=mean;
433 vNoiseRMS[i]=rms;
434
435 if (cs[i]>verySmall){
436 meanSen=vNoiseMeanSenRegions[i]/cs[i];
437 rmsSen=vNoiseRMSSenRegions[i];
438 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
439 }
440 vNoiseMeanSenRegions[i]=meanSen;
441 vNoiseRMSSenRegions[i]=rmsSen;
442 }
443}
444
445//_____________________________________________________________________________________
446void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
447{
448 //
449 // Process the Pulser information
450 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
451 //
452
453 const UInt_t infoSize=4;
454 //reset counters to error number
455 vMeanTime.ResizeTo(infoSize);
456 vMeanTime.Zero();
457 //counter
458 TVectorD c(infoSize);
459 //retrieve pulser and ALTRO data
460 if (!fPulserTmean) return;
461 //
462 //get Outliers
463 AliTPCCalROC *rocOut=0x0;
464 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
465 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
466 if (!tmeanROC) continue;
467 rocOut=fPulserOutlier->GetCalROC(isec);
468 UInt_t nchannels=tmeanROC->GetNchannels();
469 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
470 if (rocOut && rocOut->GetValue(ichannel)) continue;
471 Float_t val=tmeanROC->GetValue(ichannel);
472 Int_t type=isec/18;
473 vMeanTime[type]+=val;
474 ++c[type];
475 }
476 }
477 //calculate mean
478 for (UInt_t itype=0; itype<infoSize; ++itype){
479 if (c[itype]>0) vMeanTime[itype]/=c[itype];
480 else vMeanTime[itype]=0;
481 }
482}
483//_____________________________________________________________________________________
484void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
485{
486 //
487 // Get Values from ALTRO configuration data
488 //
489 nMasked=-1;
490 if (!fALTROMasked) return;
491 nMasked=0;
492 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
493 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
494 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
495 if (rocMasked->GetValue(ichannel)) ++nMasked;
496 }
497 }
498}
499//_____________________________________________________________________________________
500void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
501{
502 //
503 // Proces Goofie values, return statistical information of the currently set goofieArray
504 // The meaning of the entries are given below
505 /*
506 1 TPC_ANODE_I_A00_STAT
507 2 TPC_DVM_CO2
508 3 TPC_DVM_DriftVelocity
509 4 TPC_DVM_FCageHV
510 5 TPC_DVM_GainFar
511 6 TPC_DVM_GainNear
512 7 TPC_DVM_N2
513 8 TPC_DVM_NumberOfSparks
514 9 TPC_DVM_PeakAreaFar
515 10 TPC_DVM_PeakAreaNear
516 11 TPC_DVM_PeakPosFar
517 12 TPC_DVM_PeakPosNear
518 13 TPC_DVM_PickupHV
519 14 TPC_DVM_Pressure
520 15 TPC_DVM_T1_Over_P
521 16 TPC_DVM_T2_Over_P
522 17 TPC_DVM_T_Over_P
523 18 TPC_DVM_TemperatureS1
524 */
525 if (!fGoofieArray){
526 Int_t nsensors=19;
527 vecEntries.ResizeTo(nsensors);
528 vecMedian.ResizeTo(nsensors);
529 vecMean.ResizeTo(nsensors);
530 vecRMS.ResizeTo(nsensors);
531 vecEntries.Zero();
532 vecMedian.Zero();
533 vecMean.Zero();
534 vecRMS.Zero();
535 return;
536 }
537 Double_t kEpsilon=0.0000000001;
538 Double_t kBig=100000000000.;
539 Int_t nsensors = fGoofieArray->NumSensors();
540 vecEntries.ResizeTo(nsensors);
541 vecMedian.ResizeTo(nsensors);
542 vecMean.ResizeTo(nsensors);
543 vecRMS.ResizeTo(nsensors);
544 TVectorF values;
545 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
546 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
547 if (gsensor && gsensor->GetGraph()){
548 Int_t npoints = gsensor->GetGraph()->GetN();
549 // filter zeroes
550 values.ResizeTo(npoints);
551 Int_t nused =0;
552 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
553 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
554 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
555 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
556 nused++;
557 }
558 }
559 //
560 vecEntries[isensor]= nused;
561 if (nused>1){
562 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
563 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
564 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
565 }
566 }
567 }
568}
7390f655 569//_____________________________________________________________________________________
570void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
571{
572 //
573 // check the variations of the pedestal data to the reference pedestal data
574 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
575 //
576 const Int_t npar=4;
577 TVectorF vThres(npar); //thresholds
578 Int_t nActive=0; //number of active channels
579
580 //reset and set thresholds
581 pedestalDeviations.ResizeTo(npar);
582 for (Int_t i=0;i<npar;++i){
583 pedestalDeviations.GetMatrixArray()[i]=0;
584 vThres.GetMatrixArray()[i]=(i+1)*.5;
585 }
586 //check all needed data is available
587 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
588 //loop over all channels
589 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
590 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
591 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
592 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
593 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
594 UInt_t nrows=mROC->GetNrows();
595 for (UInt_t irow=0;irow<nrows;++irow){
596 UInt_t npads=mROC->GetNPads(irow);
597 for (UInt_t ipad=0;ipad<npads;++ipad){
598 //don't use masked channels;
599 if (mROC ->GetValue(irow,ipad)) continue;
600 if (mRefROC->GetValue(irow,ipad)) continue;
601 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
602 for (Int_t i=0;i<npar;++i){
603 if (deviation>vThres[i])
604 ++pedestalDeviations.GetMatrixArray()[i];
605 }
606 ++nActive;
607 }//end ipad
608 }//ind irow
609 }//end isec
610 if (nActive>0){
611 for (Int_t i=0;i<npar;++i){
612 pedestalDeviations.GetMatrixArray()[i]/=nActive;
613 }
614 }
615}
616//_____________________________________________________________________________________
617void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
618{
619 //
620 // check the variations of the noise data to the reference noise data
621 // thresholds are 5, 10, 15 and 20 percent respectively.
622 //
623 const Int_t npar=4;
624 TVectorF vThres(npar); //thresholds
625 Int_t nActive=0; //number of active channels
626
627 //reset and set thresholds
628 noiseDeviations.ResizeTo(npar);
629 for (Int_t i=0;i<npar;++i){
630 noiseDeviations.GetMatrixArray()[i]=0;
631 vThres.GetMatrixArray()[i]=(i+1)*.05;
632 }
633 //check all needed data is available
634 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
635 //loop over all channels
636 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
637 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
638 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
639 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
640 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
641 UInt_t nrows=mROC->GetNrows();
642 for (UInt_t irow=0;irow<nrows;++irow){
643 UInt_t npads=mROC->GetNPads(irow);
644 for (UInt_t ipad=0;ipad<npads;++ipad){
645 //don't use masked channels;
646 if (mROC ->GetValue(irow,ipad)) continue;
647 if (mRefROC->GetValue(irow,ipad)) continue;
648 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
649 for (Int_t i=0;i<npar;++i){
650 if (deviation>vThres[i])
651 ++noiseDeviations.GetMatrixArray()[i];
652 }
653 ++nActive;
654 }//end ipad
655 }//ind irow
656 }//end isec
657 if (nActive>0){
658 for (Int_t i=0;i<npar;++i){
659 noiseDeviations.GetMatrixArray()[i]/=nActive;
660 }
661 }
662}
663//_____________________________________________________________________________________
664void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
665 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
666{
667 //
668 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
669 // thresholds are .5, 1, 5 and 10 percent respectively.
670 //
671 //
672 const Int_t npar=4;
673 TVectorF vThres(npar); //thresholds
674 Int_t nActive=0; //number of active channels
675
676 //reset and set thresholds
677 pulserQdeviations.ResizeTo(npar);
678 for (Int_t i=0;i<npar;++i){
679 pulserQdeviations.GetMatrixArray()[i]=0;
680 }
681 npadsOutOneTB=0;
682 npadsOffAdd=0;
683 varQMean=0;
684 vThres.GetMatrixArray()[0]=.005;
685 vThres.GetMatrixArray()[1]=.01;
686 vThres.GetMatrixArray()[2]=.05;
687 vThres.GetMatrixArray()[3]=.1;
688 //check all needed data is available
689 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
690 //
691 UpdateRefPulserOutlierMap();
692 //loop over all channels
732e90a8 693 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 694 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
695 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
696 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
732e90a8 697// AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
7390f655 698 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
699 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
700 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
701 Float_t pt_mean=ptROC->GetMean(oROC);
702 UInt_t nrows=mROC->GetNrows();
703 for (UInt_t irow=0;irow<nrows;++irow){
704 UInt_t npads=mROC->GetNPads(irow);
705 for (UInt_t ipad=0;ipad<npads;++ipad){
706 //don't use masked channels;
707 if (mROC ->GetValue(irow,ipad)) continue;
708 if (mRefROC->GetValue(irow,ipad)) continue;
709 //don't user edge pads
710 if (ipad==0||ipad==npads-1) continue;
711 //data
712 Float_t pq=pqROC->GetValue(irow,ipad);
713 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
714 Float_t pt=ptROC->GetValue(irow,ipad);
715// Float_t ptRef=ptRefROC->GetValue(irow,ipad);
716 //comparisons q
717 Float_t deviation=TMath::Abs(pq/pqRef-1);
718 for (Int_t i=0;i<npar;++i){
719 if (deviation>vThres[i])
720 ++pulserQdeviations.GetMatrixArray()[i];
721 }
722 if (pqRef>11&&pq<11) ++npadsOffAdd;
723 varQMean+=pq-pqRef;
724 //comparisons t
725 if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB;
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 pulserQdeviations.GetMatrixArray()[i]/=nActive;
733 varQMean/=nActive;
734 }
735 }
736}
892226be 737//_____________________________________________________________________________________
738void AliTPCcalibDButil::UpdatePulserOutlierMap()
7390f655 739{
740 //
741 //
742 //
743 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
744}
745//_____________________________________________________________________________________
746void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
747{
748 //
749 //
750 //
751 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
752}
753//_____________________________________________________________________________________
754void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
892226be 755{
756 //
757 // Create a map that contains outliers from the Pulser calibration data.
758 // The outliers include masked channels, edge pads and pads with
759 // too large timing and charge variations.
7390f655 760 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
892226be 761 // those do not contain masked and edge pads
762 //
7390f655 763 if (!pulT||!pulQ) {
892226be 764 //reset map
7390f655 765 pulOut->Multiply(0.);
892226be 766 fNpulserOutliers=-1;
767 return;
768 }
769 AliTPCCalROC *rocMasked=0x0;
770 fNpulserOutliers=0;
771
772 //Create Outlier Map
773 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 774 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
775 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
776 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
892226be 777 if (!tmeanROC||!qmeanROC) {
778 //reset outliers in this ROC
779 outROC->Multiply(0.);
780 continue;
781 }
782 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
783// Double_t dummy=0;
784// Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
785// Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
786 UInt_t nrows=tmeanROC->GetNrows();
787 for (UInt_t irow=0;irow<nrows;++irow){
788 UInt_t npads=tmeanROC->GetNPads(irow);
789 for (UInt_t ipad=0;ipad<npads;++ipad){
790 Int_t outlier=0,masked=0;
791 Float_t q=qmeanROC->GetValue(irow,ipad);
792 Float_t t=tmeanROC->GetValue(irow,ipad);
793 //masked channels are outliers
794 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
795 //edge pads are outliers
796 if (ipad==0||ipad==npads-1) masked=1;
797 //channels with too large charge or timing deviation from the meadian are outliers
798// if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
799 if (q<fPulQminLimit && !masked) outlier=1;
800 //check for nan
801 if ( !(q<10000000) || !(t<10000000)) outlier=1;
802 outROC->SetValue(irow,ipad,outlier+masked);
803 fNpulserOutliers+=outlier;
804 }
805 }
806 }
807}
808//_____________________________________________________________________________________
2cb269df 809AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
892226be 810{
811 //
812 // Create pad time0 object from pulser and/or CE data, depending on the selected model
813 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
814 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
732e90a8 815 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
892226be 816 //
2cb269df 817 // In case model 2 is invoked - gy arival time gradient is also returned
818 //
819 gyA=0;
820 gyC=0;
892226be 821 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
822 // decide between different models
823 if (model==0||model==1){
824 TVectorD vMean;
825 if (model==1) ProcessPulser(vMean);
826 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
827 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
828 if (!rocPulTmean) continue;
829 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
830 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
831 Float_t mean=rocPulTmean->GetMean(rocOut);
832 //treat case where a whole partition is masked
833 if (mean==0) mean=rocPulTmean->GetMean();
834 if (model==1) {
835 Int_t type=isec/18;
836 mean=vMean[type];
837 }
838 UInt_t nrows=rocTime0->GetNrows();
839 for (UInt_t irow=0;irow<nrows;++irow){
840 UInt_t npads=rocTime0->GetNPads(irow);
841 for (UInt_t ipad=0;ipad<npads;++ipad){
842 Float_t time=rocPulTmean->GetValue(irow,ipad);
843 //in case of an outlier pad use the mean of the altro values.
844 //This should be the most precise guess in that case.
845 if (rocOut->GetValue(irow,ipad)) {
846 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
847 if (time==0) time=mean;
848 }
849 Float_t val=time-mean;
850 rocTime0->SetValue(irow,ipad,val);
851 }
852 }
853 }
2cb269df 854 } else if (model==2){
855 Double_t pgya,pgyc,pchi2a,pchi2c;
856 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
857 fCETmean->Add(padPulser,-1.);
732e90a8 858 TVectorD vA,vC;
859 AliTPCCalPad outCE("outCE","outCE");
860 Int_t nOut;
2cb269df 861 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
862 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
732e90a8 863// AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
2cb269df 864 if (!padFit) { delete padPulser; return 0;}
865 gyA=vA[2];
866 gyC=vC[2];
867 fCETmean->Add(padPulser,1.);
732e90a8 868 padTime0->Add(fCETmean);
2cb269df 869 padTime0->Add(padFit,-1);
870 delete padPulser;
732e90a8 871 TVectorD vFitROC;
872 TMatrixD mFitROC;
873 Float_t chi2;
874 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
875 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
876 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
877 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
878 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
879 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
880 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
881 Float_t mean=rocPulTmean->GetMean(rocOutPul);
882 if (mean==0) mean=rocPulTmean->GetMean();
883 UInt_t nrows=rocTime0->GetNrows();
884 for (UInt_t irow=0;irow<nrows;++irow){
885 UInt_t npads=rocTime0->GetNPads(irow);
886 for (UInt_t ipad=0;ipad<npads;++ipad){
887 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
888 if (rocOutCE->GetValue(irow,ipad)){
889 Float_t valOut=rocCEfit->GetValue(irow,ipad);
890 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
891 rocTime0->SetValue(irow,ipad,valOut);
892 }
893 }
894 }
895 delete rocCEfit;
896 }
897 delete padFit;
892226be 898 }
2cb269df 899 Double_t median = padTime0->GetMedian();
900 padTime0->Add(-median); // normalize to median
892226be 901 return padTime0;
902}
903//_____________________________________________________________________________________
904Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
905{
906 if (roc==0) return 0.;
907 const Int_t sector=roc->GetSector();
908 AliTPCROC *tpcRoc=AliTPCROC::Instance();
909 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
910 Float_t mean=0;
911 Int_t n=0;
912
913 //loop over a small range around the requested pad (+-10 rows/pads)
914 for (Int_t irow=row-10;irow<row+10;++irow){
915 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
916 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
917 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
918 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
919 if (altroRoc!=altroCurr) continue;
920 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
921 Float_t val=roc->GetValue(irow,ipad);
922 mean+=val;
923 ++n;
924 }
925 }
926 if (n>0) mean/=n;
927 return mean;
928}
7390f655 929//_____________________________________________________________________________________
930void AliTPCcalibDButil::SetRefFile(const char* filename)
931{
932 //
933 // load cal pad objects form the reference file
934 //
935 TDirectory *currDir=gDirectory;
936 TFile f(filename);
937 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
938 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
939 //pulser data
940 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
941 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
942 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
943 //CE data
944 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
945 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
946 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
947 //Altro data
948// fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
949// fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
950// fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
951// fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
952 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
953 f.Close();
954 currDir->cd();
955}
956
2cb269df 957
958
959
960AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad *ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
961 //
962 // Author: marian.ivanov@cern.ch
963 //
964 // Create outlier map for CE study
965 // Parameters:
966 // Return value - outlyer map
967 // noutlyersCE - number of outlyers
968 // minSignal - minimal total Q signal
969 // cutRMSMin - minimal width of the signal in respect to the median
970 // cutRMSMax - maximal width of the signal in respect to the median
971 // cutMaxDistT - maximal deviation from time median per chamber
972 //
973 // Outlyers criteria:
974 // 0. Exclude masked pads
975 // 1. Exclude first two rows in IROC and last two rows in OROC
976 // 2. Exclude edge pads
977 // 3. Exclude channels with too large variations
978 // 4. Exclude pads with too small signal
979 // 5. Exclude signal with outlyers RMS
980 // 6. Exclude channels to far from the chamber median
981 noutliersCE=0;
982 //create outlier map
983 AliTPCCalPad *out=ceOut;
984 if (!out) out= new AliTPCCalPad("outCE","outCE");
985 AliTPCCalROC *rocMasked=0x0;
986 if (!fCETmean) return 0;
987 if (!fCETrms) return 0;
988 if (!fCEQmean) return 0;
989 //
990 //loop over all channels
991 //
992 Double_t rmsMedian = fCETrms->GetMedian();
993 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
994 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
995 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
996 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
997 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
998 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
999 Double_t trocMedian = rocData->GetMedian();
1000 //
1001 if (!rocData) {
1002 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1003 rocOut->Add(1.);
1004 continue;
1005 }
1006 //
1007 //select outliers
1008 UInt_t nrows=rocData->GetNrows();
1009 for (UInt_t irow=0;irow<nrows;++irow){
1010 UInt_t npads=rocData->GetNPads(irow);
1011 for (UInt_t ipad=0;ipad<npads;++ipad){
1012 rocOut->SetValue(irow,ipad,0);
1013 Float_t valTmean=rocData->GetValue(irow,ipad);
1014 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1015 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1016 //0. exclude masked pads
1017 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1018 rocOut->SetValue(irow,ipad,1);
1019 continue;
1020 }
1021 //1. exclude first two rows in IROC and last two rows in OROC
1022 if (iroc<36){
1023 if (irow<2) rocOut->SetValue(irow,ipad,1);
1024 } else {
1025 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1026 }
1027 //2. exclude edge pads
1028 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1029 //exclude values that are exactly 0
1030 if (valTmean==0) {
1031 rocOut->SetValue(irow,ipad,1);
1032 ++noutliersCE;
1033 }
1034 //3. exclude channels with too large variations
1035 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1036 rocOut->SetValue(irow,ipad,1);
1037 ++noutliersCE;
1038 }
1039 //
1040 //4. exclude channels with too small signal
1041 if (valQmean<minSignal) {
1042 rocOut->SetValue(irow,ipad,1);
1043 ++noutliersCE;
1044 }
1045 //
1046 //5. exclude channels with too small rms
1047 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1048 rocOut->SetValue(irow,ipad,1);
1049 ++noutliersCE;
1050 }
1051 //
1052 //6. exclude channels to far from the chamber median
1053 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1054 rocOut->SetValue(irow,ipad,1);
1055 ++noutliersCE;
1056 }
1057 }
1058 }
1059 }
1060 //
1061 return out;
1062}
1063
1064
1065AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad *pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1066 //
1067 // Author: marian.ivanov@cern.ch
1068 //
1069 // Create outlier map for Pulser
1070 // Parameters:
1071 // Return value - outlyer map
1072 // noutlyersPulser - number of outlyers
1073 // cutTime - absolute cut - distance to the median of chamber
1074 // cutnRMSQ - nsigma cut from median q distribution per chamber
1075 // cutnRMSrms - nsigma cut from median rms distribution
1076 // Outlyers criteria:
1077 // 0. Exclude masked pads
1078 // 1. Exclude time outlyers (default 3 time bins)
1079 // 2. Exclude q outlyers (default 5 sigma)
1080 // 3. Exclude rms outlyers (default 5 sigma)
1081 noutliersPulser=0;
1082 AliTPCCalPad *out=pulserOut;
1083 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1084 AliTPCCalROC *rocMasked=0x0;
1085 if (!fPulserTmean) return 0;
1086 if (!fPulserTrms) return 0;
1087 if (!fPulserQmean) return 0;
1088 //
1089 //loop over all channels
1090 //
1091 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1092 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1093 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1094 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1095 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1096 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1097 //
1098 Double_t rocMedianT = rocData->GetMedian();
1099 Double_t rocMedianQ = rocPulserQ->GetMedian();
1100 Double_t rocRMSQ = rocPulserQ->GetRMS();
1101 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1102 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1103 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1104 rocOut->SetValue(ichannel,0);
1105 Float_t valTmean=rocData->GetValue(ichannel);
1106 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1107 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1108 Int_t isOut=0;
1109 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1110 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1111 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1112 rocOut->SetValue(ichannel,isOut);
1113 if (isOut) noutliersPulser++;
1114 }
1115 }
1116 return out;
1117}
1118
1119
1120AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1121 //
1122 // Author : Marian Ivanov
1123 // Create pad time0 correction map using information from the CE and from pulser
1124 //
1125 //
1126 // Return PadTime0 to be used for time0 relative alignment
1127 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1128 // using $ALICE_ROOT/TPC/script/gui application
1129 //
1130 // fitResultsA - fitParameters A side
1131 // fitResultsC - fitParameters C side
1132 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1133 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1134 //
1135 //
1136 // Algorithm:
1137 // 1. Find outlier map for CE
1138 // 2. Find outlier map for Pulser
1139 // 3. Replace outlier by median at given sector (median without outliers)
1140 // 4. Substract from the CE data pulser
1141 // 5. Fit the CE with formula
1142 // 5.1) (IROC-OROC) offset
1143 // 5.2) gx
1144 // 5.3) gy
1145 // 5.4) (lx-xmid)
1146 // 5.5) (IROC-OROC)*(lx-xmid)
1147 // 5.6) (ly/lx)^2
1148 // 6. Substract gy fit dependence from the CE data
1149 // 7. Add pulser back to CE data
1150 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1151 // 9. return CE data
1152 //
1153 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1154 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1155
1156 // fit formula
1157 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)";
1158 // output for fit formula
1159 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)";
1160 // gy part of formula
1161 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)";
1162 //
1163 //
1164 if (!fCETmean) return 0;
1165 Double_t pgya,pgyc,pchi2a,pchi2c;
1166 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1167 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1168
1169 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1170 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1171 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1172 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1173 padPulser->SetName("padPulser");
1174 padPulserOut->SetName("padPulserOut");
1175 padCE->SetName("padCE");
1176 padCEIn->SetName("padCEIn");
1177 padCEOut->SetName("padCEOut");
1178 padOut->SetName("padOut");
1179
1180 //
1181 // make combined outlyers map
1182 // and replace outlyers in maps with median for chamber
1183 //
1184 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1185 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1186 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1187 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1188 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1189 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1190 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1191 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1192 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1193 if (rocPulserOut->GetValue(ichannel)>0) {
1194 rocPulser->SetValue(ichannel,pulserMedian);
1195 rocOut->SetValue(ichannel,1);
1196 }
1197 if (rocCEOut->GetValue(ichannel)>0) {
1198 rocCE->SetValue(ichannel,ceMedian);
1199 rocOut->SetValue(ichannel,1);
1200 }
1201 }
1202 }
1203 //
1204 // remove pulser time 0
1205 //
1206 padCE->Add(padPulser,-1);
1207 //
1208 // Make fits
1209 //
1210 TMatrixD dummy;
1211 Float_t chi2Af,chi2Cf;
1212 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1213 chi2A=chi2Af;
1214 chi2C=chi2Cf;
1215 //
1216 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1217 padCEFitGY->SetName("padCEFitGy");
1218 //
1219 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1220 padCEFit->SetName("padCEFit");
1221 //
1222 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1223 padCEDiff->SetName("padCEDiff");
1224 padCEDiff->Add(padCEFit,-1.);
1225 //
1226 //
1227 padCE->Add(padCEFitGY,-1.);
1228
1229 padCE->Add(padPulser,1.);
1230 Double_t padmedian = padCE->GetMedian();
1231 padCE->Add(-padmedian); // normalize to median
1232 //
1233 // Replace outliers by fit value - median of diff per given chamber -GY fit
1234 //
1235 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1236 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1237 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1238 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1239 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1240 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1241 //
1242 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1243 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1244 if (rocOut->GetValue(ichannel)==0) continue;
1245 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1246 rocCE->SetValue(ichannel,value);
1247 }
1248 }
1249 //
1250 //
1251 if (dumpfile){
1252 //dump to the file - result can be visualized
1253 AliTPCPreprocessorOnline preprocesor;
1254 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1255 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1256 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1257 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1258 //
1259 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1260 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1261 //
1262 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1263 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1264 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1265 preprocesor.DumpToFile(dumpfile);
1266 }
1267 delete padPulser;
1268 delete padPulserOut;
1269 delete padCEIn;
1270 delete padCEOut;
1271 delete padOut;
1272 delete padCEDiff;
1273 delete padCEFitGY;
1274 return padCE;
1275}
1276
817766d5 1277
1278
1279
1280
1281Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1282 //
1283 // find the closest point to xref in x direction
1284 // return dx and value
1285 Int_t index=0;
1286 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1287 if (index<0) index=0;
1288 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1289 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1290 dx = xref-graph->GetX()[index];
1291 y = graph->GetY()[index];
1292 return index;
1293}
1294
1295
1296Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1297 //
1298 // Get the correction of the trigger offset
1299 // combining information from the laser track calibration
1300 // and from cosmic calibration
1301 //
1302 // run - run number
1303 // timeStamp - tim stamp in seconds
1304 // deltaT - integration period to calculate offset
1305 // deltaTLaser -max validity of laser data
1306 // valType - 0 - median, 1- mean
1307 //
1308 // Integration vaues are just recomendation - if not possible to get points
1309 // automatically increase the validity by factor 2
1310 // (recursive algorithm until one month of data taking)
1311 //
1312 //
1313 const Float_t kLaserCut=0.0005;
1314 const Int_t kMaxPeriod=3600*24*30*3; // 3 month max
1315 const Int_t kMinPoints=20;
1316 //
1317 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1318 if (!array) {
1319 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1320 }
1321 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1322 if (!array) return 0;
1323 //
1324 TGraphErrors *laserA[3]={0,0,0};
1325 TGraphErrors *laserC[3]={0,0,0};
1326 TGraphErrors *cosmicAll=0;
1327 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1328 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1329 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1330 //
1331 //
1332 if (!cosmicAll) return 0;
1333 Int_t nmeasC=cosmicAll->GetN();
1334 Float_t *tdelta = new Float_t[nmeasC];
1335 Int_t nused=0;
1336 for (Int_t i=0;i<nmeasC;i++){
1337 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1338 Float_t ccosmic=cosmicAll->GetY()[i];
1339 Double_t yA=0,yC=0,dA=0,dC=0;
1340 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1341 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1342 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1343 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1344 //
1345 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1346 Float_t claser=0;
1347 if (TMath::Abs(yA-yC)<kLaserCut) {
1348 claser=(yA-yC)*0.5;
1349 }else{
1350 if (i%2==0) claser=yA;
1351 if (i%2==1) claser=yC;
1352 }
1353 tdelta[nused]=ccosmic-claser;
1354 nused++;
1355 }
1356 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1357 Double_t median = TMath::Median(nused,tdelta);
1358 Double_t mean = TMath::Mean(nused,tdelta);
1359 delete tdelta;
1360 return (valType==0) ? median:mean;
1361}
1362
a23ba1c3 1363Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
817766d5 1364 //
1365 // Get the correction of the drift velocity
1366 // combining information from the laser track calibration
1367 // and from cosmic calibration
1368 //
a23ba1c3 1369 // dist - return value - distance to closest point in graph
817766d5 1370 // run - run number
1371 // timeStamp - tim stamp in seconds
1372 // deltaT - integration period to calculate time0 offset
1373 // deltaTLaser -max validity of laser data
1374 // valType - 0 - median, 1- mean
1375 //
1376 // Integration vaues are just recomendation - if not possible to get points
1377 // automatically increase the validity by factor 2
1378 // (recursive algorithm until one month of data taking)
1379 //
1380 //
1381 //
1382 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1383 if (!array) {
1384 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1385 }
1386 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1387 if (!array) return 0;
1388 TGraphErrors *cosmicAll=0;
1389 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1390 if (!cosmicAll) return 0;
a23ba1c3 1391 Double_t grY=0;
1392 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1393
817766d5 1394 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
a23ba1c3 1395 Double_t vcosmic= AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1e722a63 1396 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
817766d5 1397 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1398 return vcosmic+t0;
1399
1400 /*
1401 Example usage:
1402
1403 Int_t run=89000
1404 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1405 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1406 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1407 //
1408 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1409 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1410 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1411 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1412
1413 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1414 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1415
1416 */
1417
1418}
1419
a23ba1c3 1420Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1421 //
1422 // Get the correction of the drift velocity using the laser tracks calbration
1423 //
1424 // run - run number
1425 // timeStamp - tim stamp in seconds
1426 // deltaT - integration period to calculate time0 offset
1427 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1428 // Note in case no data form both A and C side - the value from active side used
1429 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1430 TGraphErrors *grlaserA=0;
1431 TGraphErrors *grlaserC=0;
1432 Double_t vlaserA=0, vlaserC=0;
1433 if (!array) return 0;
1434 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1435 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1436 Double_t deltaY;
1437 if (grlaserA) {
1438 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1439 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1440 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1441 }
1442 if (grlaserC) {
1443 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1444 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1445 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1446 }
1447 if (side==0) return vlaserA;
1448 if (side==1) return vlaserC;
1449 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1450 if (!grlaserA) return vlaserC;
1451 if (!grlaserC) return vlaserA;
1452 return mdrift;
1453}
1454
1455
1456Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1457 //
1458 // Get the correction of the drift velocity using the CE laser data
1459 // combining information from the CE, laser track calibration
1460 // and P/T calibration
1461 //
1462 // run - run number
1463 // timeStamp - tim stamp in seconds
1464 // deltaT - integration period to calculate time0 offset
1465 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1466 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
1467 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1468 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1469 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1470 //
1471 //
1472 Double_t corrPTA = 0, corrPTC=0;
1473 Double_t ltime0A = 0, ltime0C=0;
1474 Double_t gry=0;
1475 Double_t corrA=0, corrC=0;
1476 Double_t timeA=0, timeC=0;
1477 TGraph *graphA = (TGraph*)arrT->At(72);
1478 TGraph *graphC = (TGraph*)arrT->At(73);
1479 if (!graphA && !graphC) return 0.;
1480 if (graphA &&graphA->GetN()>0) {
1481 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1482 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1483 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1484 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1485 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 1486 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 1487 corrA-=corrPTA;
1488 }
1489 if (graphC&&graphC->GetN()>0){
1490 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1491 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1492 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1493 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1494 if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 1495 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 1496 corrC-=corrPTC;
1497 }
1498
1499 if (side ==0 ) return corrA;
1500 if (side ==1 ) return corrC;
1501 Double_t corrM= (corrA+corrC)*0.5;
1502 if (!graphA) corrM=corrC;
1503 if (!graphC) corrM=corrA;
1504 return corrM;
1505}
1506
1507
1508
1509
1510Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1511 //
1512 // VERY obscure method - we need something in framework
1513 // Find the TPC runs with temperature OCDB entry
1514 // cache the start and end of the run
1515 //
1516 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1517 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1518 if (!storage) return 0;
1519 TString path=storage->GetURI();
1520 TString runsT;
1521 {
1522 TString command;
1523 if (path.Contains("local")){ // find the list if local system
1524 path.ReplaceAll("local://","");
1525 path+="TPC/Calib/Temperature";
1526 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
1527 }
1528 runsT=gSystem->GetFromPipe(command);
1529 }
1530 TObjArray *arr= runsT.Tokenize("r");
1531 if (!arr) return 0;
1532 //
1533 TArrayI indexes(arr->GetEntries());
1534 TArrayI runs(arr->GetEntries());
1535 Int_t naccept=0;
1536 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1537 Int_t irunN = atoi(arr->At(irun)->GetName());
1538 if (irunN<startRun) continue;
1539 if (irunN>stopRun) continue;
1540 runs[naccept]=irunN;
1541 naccept++;
1542 }}
1543 fRuns.Set(naccept);
1544 fRunsStart.Set(fRuns.fN);
1545 fRunsStop.Set(fRuns.fN);
1546 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1547 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
1548
1549 //
1550 AliCDBEntry * entry = 0;
1551 {for (Int_t irun=0;irun<fRuns.fN; irun++){
1552 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
1553 if (!entry) continue;
1554 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
1555 if (!tmpRun) continue;
1556 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
1557 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
1558 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
1559 }}
1560 return fRuns.fN;
1561}
1562
1563
1564Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
1565 //
1566 // binary search - find the run for given time stamp
1567 //
1568 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
1569 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
1570 Int_t cindex = -1;
1571 for (Int_t index=index0; index<=index1; index++){
1572 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
1573 if (debug) {
1574 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
1575 }
1576 }
1577 if (cindex<0) cindex =(index0+index1)/2;
1578 if (cindex<0) {
1579 return 0;
1580 }
1581 return fRuns[cindex];
1582}
1583
1584
1585
1586
1587
1588TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
1589 //
1590 // filter outlyer measurement
1591 // Only points around median +- sigmaCut filtered
1592 if (!graph) return 0;
1593 Int_t kMinPoints=2;
1594 Int_t npoints0 = graph->GetN();
1595 Int_t npoints=0;
1596 Float_t rmsY=0;
1597 Double_t *outx=new Double_t[npoints0];
1598 Double_t *outy=new Double_t[npoints0];
1599 //
1600 //
1601 if (npoints0<kMinPoints) return 0;
1602 for (Int_t iter=0; iter<3; iter++){
1603 npoints=0;
1604 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1605 if (graph->GetY()[ipoint]==0) continue;
1606 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
1607 outx[npoints] = graph->GetX()[ipoint];
1608 outy[npoints] = graph->GetY()[ipoint];
1609 npoints++;
1610 }
1611 if (npoints<=1) break;
1612 medianY =TMath::Median(npoints,outy);
1613 rmsY =TMath::RMS(npoints,outy);
1614 }
1615 TGraph *graphOut=0;
1616 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1617 return graphOut;
1618}
1619
1620
1621TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
1622 //
1623 // filter outlyer measurement
1624 // Only points around median +- cut filtered
1625 if (!graph) return 0;
1626 Int_t kMinPoints=2;
1627 Int_t npoints0 = graph->GetN();
1628 Int_t npoints=0;
1629 Float_t rmsY=0;
1630 Double_t *outx=new Double_t[npoints0];
1631 Double_t *outy=new Double_t[npoints0];
1632 //
1633 //
1634 if (npoints0<kMinPoints) return 0;
1635 for (Int_t iter=0; iter<3; iter++){
1636 npoints=0;
1637 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1638 if (graph->GetY()[ipoint]==0) continue;
1639 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
1640 outx[npoints] = graph->GetX()[ipoint];
1641 outy[npoints] = graph->GetY()[ipoint];
1642 npoints++;
1643 }
1644 if (npoints<=1) break;
1645 medianY =TMath::Median(npoints,outy);
1646 rmsY =TMath::RMS(npoints,outy);
1647 }
1648 TGraph *graphOut=0;
1649 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1650 return graphOut;
1651}
1652
1653
1654
1655TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY){
1656 //
1657 // filter outlyer measurement
1658 // Only points with normalized errors median +- sigmaCut filtered
1659 //
1660 Int_t kMinPoints=10;
1661 Int_t npoints0 = graph->GetN();
1662 Int_t npoints=0;
1663 Float_t medianErr=0, rmsErr=0;
1664 Double_t *outx=new Double_t[npoints0];
1665 Double_t *outy=new Double_t[npoints0];
1666 Double_t *erry=new Double_t[npoints0];
1667 Double_t *nerry=new Double_t[npoints0];
1668 Double_t *errx=new Double_t[npoints0];
1669 //
1670 //
1671 if (npoints0<kMinPoints) return 0;
1672 for (Int_t iter=0; iter<3; iter++){
1673 npoints=0;
1674 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1675 nerry[npoints] = graph->GetErrorY(ipoint);
1676 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
1677 erry[npoints] = graph->GetErrorY(ipoint);
1678 outx[npoints] = graph->GetX()[ipoint];
1679 outy[npoints] = graph->GetY()[ipoint];
1680 errx[npoints] = graph->GetErrorY(ipoint);
1681 npoints++;
1682 }
1683 if (npoints==0) break;
1684 medianErr=TMath::Median(npoints,erry);
1685 medianY =TMath::Median(npoints,outy);
1686 rmsErr =TMath::RMS(npoints,erry);
1687 }
1688 TGraphErrors *graphOut=0;
1689 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
1690 delete []outx;
1691 delete []outy;
1692 delete []errx;
1693 delete []erry;
1694 return graphOut;
1695}
1696
1697
1698void AliTPCcalibDButil::Sort(TGraph *graph){
1699 //
1700 // sort array - neccessay for approx
1701 //
1702 Int_t npoints = graph->GetN();
1703 Int_t *indexes=new Int_t[npoints];
1704 Double_t *outx=new Double_t[npoints];
1705 Double_t *outy=new Double_t[npoints];
1706 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
1707 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
1708 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
1709 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
1710 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
1711
1712}
1713void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
1714 //
1715 // smmoth graph - mean on the interval
1716 //
1717 Sort(graph);
1718 Int_t npoints = graph->GetN();
1719 Double_t *outy=new Double_t[npoints];
1720 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1721 Double_t lx=graph->GetX()[ipoint];
1722 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
1723 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
1724 if (index0<0) index0=0;
1725 if (index1>=npoints-1) index1=npoints-1;
1726 if ((index1-index0)>1){
1727 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
1728 }else{
1729 outy[ipoint]=graph->GetY()[ipoint];
1730 }
1731 }
1732 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1733 graph->GetY()[ipoint] = outy[ipoint];
1734 }
1735 delete[] outy;
1736}
1737
1738Double_t AliTPCcalibDButil::EvalGraphConst(TGraph *graph, Double_t xref){
1739 //
1740 // Use constant interpolation outside of range
1741 //
1742 if (!graph) {
1743 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
1744 return 0;
1745 }
1746 if (graph->GetN()<1){
1747 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
1748 return 0;
1749 }
1750 if (xref<graph->GetX()[0]) return graph->GetY()[0];
1751 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
1752 return graph->Eval( xref);
1753}
1754
1e722a63 1755Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
1756 //
1757 // Filter DCS sensor information
1758 // ymin - minimal value
1759 // ymax - max value
1760 // maxdy - maximal deirivative
1761 // sigmaCut - cut on values and derivative in terms of RMS distribution
1762 // Return value - accepted fraction
1763 //
1764 // Algorithm:
1765 //
1766 // 0. Calculate median and rms of values in specified range
1767 // 1. Filter out outliers - median+-sigmaCut*rms
1768 // values replaced by median
1769 //
1770 AliSplineFit * fit = sensor->GetFit();
1771 if (!fit) return 0.;
1772 Int_t nknots = fit->GetKnots();
1773 if (nknots==0) {
1774 delete fit;
1775 sensor->SetFit(0);
1776 return 0;
1777 }
1778 //
1779 Double_t *yin0 = new Double_t[nknots];
1780 Double_t *yin1 = new Double_t[nknots];
1781 Int_t naccept=0;
1782
1783 for (Int_t iknot=0; iknot< nknots; iknot++){
1784 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
1785 yin0[naccept] = fit->GetY0()[iknot];
1786 yin1[naccept] = fit->GetY1()[iknot];
1787 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
1788 naccept++;
1789 }
1790 }
1791 if (naccept<1) {
1792 delete fit;
1793 sensor->SetFit(0);
1794 return 0.;
1795 }
1796 Double_t medianY0=0, medianY1=0;
1797 Double_t rmsY0 =0, rmsY1=0;
1798 medianY0 = TMath::Median(naccept, yin0);
1799 medianY1 = TMath::Median(naccept, yin1);
1800 rmsY0 = TMath::RMS(naccept, yin0);
1801 rmsY1 = TMath::RMS(naccept, yin1);
1802 naccept=0;
1803 //
1804 // 1. Filter out outliers - median+-sigmaCut*rms
1805 // values replaced by median
1806 // if replaced the derivative set to 0
1807 //
1808 for (Int_t iknot=0; iknot< nknots; iknot++){
1809 Bool_t isOK=kTRUE;
1810 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
1811 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
1812 if (nknots<2) fit->GetY1()[iknot]=0;
1813 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
1814 if (!isOK){
1815 fit->GetY0()[iknot]=medianY0;
1816 fit->GetY1()[iknot]=0;
1817 }else{
1818 naccept++;
1819 }
1820 }
1821 delete [] yin0;
1822 delete [] yin1;
1823 return Float_t(naccept)/Float_t(nknots);
1824}
1825
1826Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
1827 //
1828 // Filter temperature array
1829 // tempArray - array of temperatures -
1830 // ymin - minimal accepted temperature - default 15
1831 // ymax - maximal accepted temperature - default 22
1832 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
1833 // return value - fraction of filtered sensors
1834 const Double_t kMaxDy=0.1;
1835 Int_t nsensors=tempArray->NumSensors();
1836 if (nsensors==0) return 0.;
1837 Int_t naccept=0;
1838 for (Int_t isensor=0; isensor<nsensors; isensor++){
1839 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
1840 if (!sensor) continue;
1841 //printf("%d\n",isensor);
1842 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
1843 if (sensor->GetFit()==0){
1844 delete sensor;
1845 tempArray->RemoveSensorNum(isensor);
1846 }else{
1847 naccept++;
1848 }
1849 }
1850 return Float_t(naccept)/Float_t(nsensors);
1851}
1852
a23ba1c3 1853
1854void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector *pcstream){
1855 //
1856 // Filter CE data
1e722a63 1857 // Input parameters:
1858 // deltaT - smoothing window (in seconds)
1859 // cutAbs - max distance of the time info to the median (in time bins)
1860 // cutSigma - max distance (in the RMS)
1861 // pcstream - optional debug streamer to store original and filtered info
1862 // Hardwired parameters:
1863 // kMinPoints =10; // minimal number of points to define the CE
1864 // kMinSectors=12; // minimal number of sectors to define sideCE
1865 // Algorithm:
1866 // 0. Filter almost emty graphs (kMinPoints=10)
1867 // 1. calculate median and RMS per side
1868 // 2. Filter graphs - in respect with side medians
1869 // - cutAbs and cutDelta used
1870 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
1871 // 4. Calculate mean for A side and C side
1872 //
1873 const Int_t kMinPoints =10; // minimal number of points to define the CE
1874 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
1875 const Int_t kMinTime =400; // minimal arrival time of CE
1876 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
a23ba1c3 1877 Double_t medianY=0;
1e722a63 1878 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
a23ba1c3 1879 if (!cearray) return;
1e722a63 1880 Double_t tmin=-1;
1881 Double_t tmax=-1;
1882 //
1883 //
a23ba1c3 1884 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
1885 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernPressure");
1886 if ( tempMapCE && cavernPressureCE){
1e722a63 1887 //
1888 Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
1889 FilterSensor(cavernPressureCE,960,1050,10, 5.);
1890 if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
1891 if (isOK) {
1892 // recalculate P/T correction map for time of the CE
1893 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
1894 driftCalib->SetName("driftPTCE");
1895 driftCalib->SetTitle("driftPTCE");
1896 cearray->AddLast(driftCalib);
1897 }
a23ba1c3 1898 }
1e722a63 1899 //
1900 // 0. Filter almost emty graphs
1901 //
a23ba1c3 1902
1e722a63 1903 for (Int_t i=0; i<72;i++){
a23ba1c3 1904 TGraph *graph= (TGraph*)arrT->At(i);
1905 if (!graph) continue;
1906 if (graph->GetN()<kMinPoints){
1907 arrT->AddAt(0,i);
1908 delete graph; // delete empty graph
1909 continue;
1910 }
1e722a63 1911 if (tmin<0) tmin = graph->GetX()[0];
1912 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
1913 //
1914 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
1915 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
1916 }
1917 //
1918 // 1. calculate median and RMS per side
1919 //
1920 TArrayF arrA(100000), arrC(100000);
1921 Int_t nA=0, nC=0;
1922 Double_t medianA=0, medianC=0;
1923 Double_t rmsA=0, rmsC=0;
1924 for (Int_t isec=0; isec<72;isec++){
1925 TGraph *graph= (TGraph*)arrT->At(isec);
1926 if (!graph) continue;
1927 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1928 if (graph->GetY()[ipoint]<kMinTime) continue;
1929 if (nA>=arrA.fN) arrA.Set(nA*2);
1930 if (nC>=arrC.fN) arrC.Set(nC*2);
1931 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
1932 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
1933 }
1934 }
1935 if (nA>0){
1936 medianA=TMath::Median(nA,arrA.fArray);
1937 rmsA =TMath::RMS(nA,arrA.fArray);
1938 }
1939 if (nC>0){
1940 medianC=TMath::Median(nC,arrC.fArray);
1941 rmsC =TMath::RMS(nC,arrC.fArray);
1942 }
1943 //
1944 // 2. Filter graphs - in respect with side medians
1945 //
1946 TArrayD vecX(100000), vecY(100000);
1947 for (Int_t isec=0; isec<72;isec++){
1948 TGraph *graph= (TGraph*)arrT->At(isec);
1949 if (!graph) continue;
1950 Double_t median = (isec%36<18) ? medianA: medianC;
1951 Double_t rms = (isec%36<18) ? rmsA: rmsC;
1952 Int_t naccept=0;
1953 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1954 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
1955 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
1956 vecX[naccept]= graph->GetX()[ipoint];
1957 vecY[naccept]= graph->GetY()[ipoint];
1958 naccept++;
1959 }
1960 if (naccept<kMinPoints){
1961 arrT->AddAt(0,isec);
1962 delete graph; // delete empty graph
1963 continue;
1964 }
1965 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
1966 delete graph;
1967 arrT->AddAt(graph2,isec);
1968 }
1969 //
1970 // 3. Cut in respect wit the graph median
1971 //
1972 for (Int_t i=0; i<72;i++){
1973 TGraph *graph= (TGraph*)arrT->At(i);
1974 if (!graph) continue;
1975 //
1976 // filter in range
1977 //
1978 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
a23ba1c3 1979 if (!graphTS0) continue;
1980 if (graphTS0->GetN()<kMinPoints) {
1981 delete graphTS0;
1982 delete graph;
1983 arrT->AddAt(0,i);
1984 continue;
1985 }
1e722a63 1986 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
a23ba1c3 1987 graphTS->Sort();
1988 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
1989 if (pcstream){
1990 Int_t run = AliTPCcalibDB::Instance()->GetRun();
1991 (*pcstream)<<"filterCE"<<
1992 "run="<<run<<
1993 "isec="<<i<<
1994 "mY="<<medianY<<
1995 "graph.="<<graph<<
1996 "graphTS0.="<<graphTS0<<
1997 "graphTS.="<<graphTS<<
1998 "\n";
1999 }
2000 delete graphTS0;
2001 if (!graphTS) continue;
2002 arrT->AddAt(graphTS,i);
2003 delete graph;
2004 }
1e722a63 2005 //
2006 // Recalculate the mean time A side C side
2007 //
2008 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2009 Int_t meanPoints=(nA+nC)/72; // mean number of points
2010 for (Int_t itime=0; itime<200; itime++){
2011 nA=0, nC=0;
2012 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2013 for (Int_t i=0; i<72;i++){
2014 TGraph *graph= (TGraph*)arrT->At(i);
2015 if (!graph) continue;
2016 if (graph->GetN()<(meanPoints/4)) continue;
2017 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2018 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2019 }
2020 xA[itime]=time;
2021 xC[itime]=time;
2022 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2023 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2024 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2025 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2026 }
2027 //
2028 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2029 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2030 if (pcstream){
2031 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2032 (*pcstream)<<"filterAC"<<
2033 "run="<<run<<
2034 "nA="<<nA<<
2035 "nC="<<nC<<
2036 "rmsTA="<<rmsTA<<
2037 "rmsTC="<<rmsTC<<
2038 "\n";
2039 }
2040 //
2041 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2042 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2043 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2044 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2045 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2046 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2047 delete grA;
2048 delete grC;
2049 if (nA<kMinSectors) arrT->AddAt(0,72);
2050 else arrT->AddAt(graphTSA,72);
2051 if (nC<kMinSectors) arrT->AddAt(0,73);
2052 else arrT->AddAt(graphTSC,73);
a23ba1c3 2053}
2054
2055
2056void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector *pcstream){
2057 //
2058 // Filter Drift velocity measurement using the tracks
2059 // 0. remove outlyers - error based
2060 // cutSigma
2061 //
2062 //
2063 const Int_t kMinPoints=1; // minimal number of points to define value
2064 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2065 Double_t medianY=0;
2066 if (!arrT) return;
2067 for (Int_t i=0; i<arrT->GetEntries();i++){
2068 TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
2069 if (!graph) continue;
2070 if (graph->GetN()<kMinPoints){
2071 delete graph;
2072 arrT->AddAt(0,i);
2073 continue;
2074 }
2075 TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2076 if (!graph2) {
2077 delete graph; arrT->AddAt(0,i); continue;
2078 }
2079 if (graph2->GetN()<1) {
2080 delete graph; arrT->AddAt(0,i); continue;
2081 }
2082 graph2->SetName(graph->GetName());
2083 graph2->SetTitle(graph->GetTitle());
2084 arrT->AddAt(graph2,i);
2085 if (pcstream){
2086 (*pcstream)<<"filterTracks"<<
2087 "run="<<run<<
2088 "isec="<<i<<
2089 "mY="<<medianY<<
2090 "graph.="<<graph<<
2091 "graph2.="<<graph2<<
2092 "\n";
2093 }
2094 delete graph;
2095 }
2096}
2097
2098
2099void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, TTreeSRedirector *pcstream){
2100 //
2101 // Filter Goofie data
2102 // 0. remove outlyers - cutSigma around median
2103 // 1. smooth the graphs - deltaT - smoothing window
2104 //
2105 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2106 //
2107 //
2108 AliDCSSensor *sensor = goofieArray->GetSensor(isensor);
2109 Double_t medianY;
2110 if (sensor && sensor->GetGraph()){
2111 TGraph * graph = sensor->GetGraph();
2112 if (isensor==3 && graph->GetN()>1){ // drift velocity
1e722a63 2113 TGraph * graphv = FilterGraphMedianAbs(graph,0.2,medianY);
a23ba1c3 2114 delete graph;
2115 graph=graphv;
2116 }
1e722a63 2117 TGraph * graph2 = FilterGraphMedian(graph,cutSigma,medianY);
a23ba1c3 2118 if (!graph2) continue;
2119 AliTPCcalibDButil::SmoothGraph(graph2,deltaT);
2120 if (pcstream){
2121 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2122 (*pcstream)<<"filterG"<<
2123 "run="<<run<<
2124 "isensor="<<isensor<<
2125 "mY="<<medianY<<
2126 "graph.="<<graph<<
2127 "graph2.="<<graph2<<
2128 "\n";
2129 }
2130 if (graph2->GetN()<2) {delete graph2; continue;}
2131 delete graph;
2132 sensor->SetGraph(graph2);
2133 }
2134 }
2135}
2136
2137
2138
2139
2140
2141
2142
2143Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2144 //
2145 //
2146 // get laser time offset
2147 // median around timeStamp+-deltaT
2148 // QA - chi2 needed for later usage - to be added
2149 // - currently cut on error
2150 //
2151 Int_t kMinPoints=1;
2152 Double_t kMinDelay=0.01;
2153 Double_t kMinDelayErr=0.0001;
2154
2155 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2156 if (!array) return 0;
2157 TGraphErrors *tlaser=0;
2158 if (array){
2159 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2160 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2161 }
2162 if (!tlaser) return 0;
2163 Int_t npoints0= tlaser->GetN();
2164 if (npoints0==0) return 0;
2165 Double_t *xlaser = new Double_t[npoints0];
2166 Double_t *ylaser = new Double_t[npoints0];
2167 Int_t npoints=0;
2168 for (Int_t i=0;i<npoints0;i++){
2169 //printf("%d\n",i);
2170 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2171 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2172 xlaser[npoints]=tlaser->GetX()[npoints];
2173 ylaser[npoints]=tlaser->GetY()[npoints];
2174 npoints++;
2175 }
2176 //
2177 //
2178 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2179 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2180 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2181 if (index0<0) index0=0;
2182 if (index1>=npoints-1) index1=npoints-1;
2183 if (index1-index0<kMinPoints) return 0;
2184 //
2185 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2186 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2187 delete [] ylaser;
2188 delete [] xlaser;
2189 return mean;
2190}