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