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