Adding AliTPCComparisonPID class
[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>
892226be 33
34#include <AliDCSSensorArray.h>
35#include <AliDCSSensor.h>
36#include "AliTPCcalibDB.h"
37#include "AliTPCCalPad.h"
38#include "AliTPCCalROC.h"
39#include "AliTPCROC.h"
40#include "AliTPCmapper.h"
41#include "AliTPCParam.h"
6e7d7dc4 42#include "AliTPCCalibRaw.h"
817766d5 43#include "TGraphErrors.h"
892226be 44
45#include "AliTPCcalibDButil.h"
2cb269df 46#include "AliTPCPreprocessorOnline.h"
892226be 47
48ClassImp(AliTPCcalibDButil)
49AliTPCcalibDButil::AliTPCcalibDButil() :
50 TObject(),
51 fCalibDB(AliTPCcalibDB::Instance()),
52 fPadNoise(0x0),
53 fPedestals(0x0),
54 fPulserTmean(0x0),
55 fPulserTrms(0x0),
56 fPulserQmean(0x0),
57 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
58 fCETmean(0x0),
59 fCETrms(0x0),
60 fCEQmean(0x0),
61 fALTROMasked(0x0),
6e7d7dc4 62 fCalibRaw(0x0),
7390f655 63 fRefPadNoise(0x0),
64 fRefPedestals(0x0),
65 fRefPulserTmean(0x0),
66 fRefPulserTrms(0x0),
67 fRefPulserQmean(0x0),
68 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
69 fRefCETmean(0x0),
70 fRefCETrms(0x0),
71 fRefCEQmean(0x0),
72 fRefALTROMasked(0x0),
73 fRefCalibRaw(0x0),
892226be 74 fGoofieArray(0x0),
75 fMapper(new AliTPCmapper(0x0)),
76 fNpulserOutliers(-1),
7390f655 77 fIrocTimeOffset(0),
892226be 78 fCETmaxLimitAbs(1.5),
79 fPulTmaxLimitAbs(1.5),
80 fPulQmaxLimitAbs(5),
81 fPulQminLimit(11)
82{
83 //
84 // Default ctor
85 //
86}
87//_____________________________________________________________________________________
88AliTPCcalibDButil::~AliTPCcalibDButil()
89{
90 //
91 // dtor
92 //
93 delete fPulserOutlier;
7390f655 94 delete fRefPulserOutlier;
892226be 95 delete fMapper;
7390f655 96 if (fRefPadNoise) delete fRefPadNoise;
97 if (fRefPedestals) delete fRefPedestals;
98 if (fRefPulserTmean) delete fRefPulserTmean;
99 if (fRefPulserTrms) delete fRefPulserTrms;
100 if (fRefPulserQmean) delete fRefPulserQmean;
101 if (fRefCETmean) delete fRefCETmean;
102 if (fRefCETrms) delete fRefCETrms;
103 if (fRefCEQmean) delete fRefCEQmean;
104 if (fRefALTROMasked) delete fRefALTROMasked;
105 if (fRefCalibRaw) delete fRefCalibRaw;
106
892226be 107}
108//_____________________________________________________________________________________
109void AliTPCcalibDButil::UpdateFromCalibDB()
110{
111 //
112 // Update pointers from calibDB
113 //
114 fPadNoise=fCalibDB->GetPadNoise();
115 fPedestals=fCalibDB->GetPedestals();
116 fPulserTmean=fCalibDB->GetPulserTmean();
117 fPulserTrms=fCalibDB->GetPulserTrms();
118 fPulserQmean=fCalibDB->GetPulserQmean();
119 fCETmean=fCalibDB->GetCETmean();
120 fCETrms=fCalibDB->GetCETrms();
121 fCEQmean=fCalibDB->GetCEQmean();
122 fALTROMasked=fCalibDB->GetALTROMasked();
123 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
6e7d7dc4 124 fCalibRaw=fCalibDB->GetCalibRaw();
892226be 125 UpdatePulserOutlierMap();
126}
127//_____________________________________________________________________________________
7390f655 128void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
2cb269df 129 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad *outCE)
892226be 130{
131 //
132 // Process the CE data for this run
133 // the return TVectorD arrays contian the results of the fit
134 // noutliersCE contains the number of pads marked as outliers,
135 // not including masked and edge pads
136 //
137
138 //retrieve CE and ALTRO data
139 if (!fCETmean){
140 TString fitString(fitFormula);
141 fitString.ReplaceAll("++","#");
142 Int_t ndim=fitString.CountChar('#')+2;
143 fitResultsA.ResizeTo(ndim);
144 fitResultsC.ResizeTo(ndim);
145 fitResultsA.Zero();
146 fitResultsC.Zero();
147 noutliersCE=-1;
148 return;
149 }
150 noutliersCE=0;
151 //create outlier map
7390f655 152 AliTPCCalPad *out=0;
153 if (outCE) out=outCE;
154 else out=new AliTPCCalPad("outCE","outCE");
892226be 155 AliTPCCalROC *rocMasked=0x0;
156 //loop over all channels
157 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
158 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
159 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
7390f655 160 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
6e7d7dc4 161 if (!rocData) {
162 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
7390f655 163 rocOut->Add(1.);
6e7d7dc4 164 continue;
165 }
892226be 166 //add time offset to IROCs
167 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
168 rocData->Add(fIrocTimeOffset);
169 //select outliers
170 UInt_t nrows=rocData->GetNrows();
171 for (UInt_t irow=0;irow<nrows;++irow){
172 UInt_t npads=rocData->GetNPads(irow);
173 for (UInt_t ipad=0;ipad<npads;++ipad){
7390f655 174 rocOut->SetValue(irow,ipad,0);
892226be 175 //exclude masked pads
176 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
7390f655 177 rocOut->SetValue(irow,ipad,1);
892226be 178 continue;
179 }
732e90a8 180 //exclude first two rows in IROC and last two rows in OROC
181 if (iroc<36){
182 if (irow<2) rocOut->SetValue(irow,ipad,1);
183 } else {
184 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
185 }
892226be 186 //exclude edge pads
7390f655 187 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
188 Float_t valTmean=rocData->GetValue(irow,ipad);
892226be 189 //exclude values that are exactly 0
190 if (valTmean==0) {
7390f655 191 rocOut->SetValue(irow,ipad,1);
892226be 192 ++noutliersCE;
193 }
194 // exclude channels with too large variations
195 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
7390f655 196 rocOut->SetValue(irow,ipad,1);
892226be 197 ++noutliersCE;
198 }
199 }
200 }
201 }
202 //perform fit
203 TMatrixD dummy;
2cb269df 204 Float_t chi2Af,chi2Cf;
205 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
206 chi2A=chi2Af;
207 chi2C=chi2Cf;
7390f655 208 if (!outCE) delete out;
892226be 209}
210//_____________________________________________________________________________________
211void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
212 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
213 Float_t &driftTimeA, Float_t &driftTimeC )
214{
215 //
216 // Calculate statistical information from the CE graphs for drift time and charge
217 //
218
219 //reset arrays
220 vecTEntries.ResizeTo(72);
221 vecTMean.ResizeTo(72);
222 vecTRMS.ResizeTo(72);
223 vecTMedian.ResizeTo(72);
224 vecQEntries.ResizeTo(72);
225 vecQMean.ResizeTo(72);
226 vecQRMS.ResizeTo(72);
227 vecQMedian.ResizeTo(72);
228 vecTEntries.Zero();
229 vecTMean.Zero();
230 vecTRMS.Zero();
231 vecTMedian.Zero();
232 vecQEntries.Zero();
233 vecQMean.Zero();
234 vecQRMS.Zero();
235 vecQMedian.Zero();
236 driftTimeA=0;
237 driftTimeC=0;
238 TObjArray *arrT=fCalibDB->GetCErocTtime();
239 TObjArray *arrQ=fCalibDB->GetCErocQtime();
240 if (arrT){
241 for (Int_t isec=0;isec<74;++isec){
242 TGraph *gr=(TGraph*)arrT->At(isec);
243 if (!gr) continue;
244 TVectorD values;
245 Int_t npoints = gr->GetN();
246 values.ResizeTo(npoints);
247 Int_t nused =0;
6e7d7dc4 248 //skip first points, theres always some problems with finding the CE position
249 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
250 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
892226be 251 values[nused]=gr->GetY()[ipoint];
252 nused++;
253 }
254 }
255 //
256 if (isec<72) vecTEntries[isec]= nused;
257 if (nused>1){
258 if (isec<72){
259 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
260 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
261 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
262 } else if (isec==72){
263 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
264 } else if (isec==73){
265 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
266 }
267 }
268 }
269 }
270 if (arrQ){
6e7d7dc4 271 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
892226be 272 TGraph *gr=(TGraph*)arrQ->At(isec);
273 if (!gr) continue;
274 TVectorD values;
275 Int_t npoints = gr->GetN();
276 values.ResizeTo(npoints);
277 Int_t nused =0;
278 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
817766d5 279 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
892226be 280 values[nused]=gr->GetY()[ipoint];
281 nused++;
282 }
283 }
284 //
6e7d7dc4 285 vecQEntries[isec]= nused;
892226be 286 if (nused>1){
287 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
288 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
289 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
290 }
291 }
292 }
293}
294
295//_____________________________________________________________________________________
296void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
297 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
298 Int_t &nonMaskedZero)
299{
300 //
301 // process noise data
302 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
303 // OROCs small pads [2] and OROCs large pads [3]
304 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
305 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
306 //
307
308 //set proper size and reset
309 const UInt_t infoSize=4;
310 vNoiseMean.ResizeTo(infoSize);
311 vNoiseMeanSenRegions.ResizeTo(infoSize);
312 vNoiseRMS.ResizeTo(infoSize);
313 vNoiseRMSSenRegions.ResizeTo(infoSize);
314 vNoiseMean.Zero();
315 vNoiseMeanSenRegions.Zero();
316 vNoiseRMS.Zero();
317 vNoiseRMSSenRegions.Zero();
318 nonMaskedZero=0;
319 //counters
320 TVectorD c(infoSize);
321 TVectorD cs(infoSize);
322 //tpc parameters
323 AliTPCParam par;
324 par.Update();
325 //retrieve noise and ALTRO data
326 if (!fPadNoise) return;
327 AliTPCCalROC *rocMasked=0x0;
328 //create IROC, OROC1, OROC2 and sensitive region masks
329 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
330 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
331 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
332 UInt_t nrows=noiseROC->GetNrows();
333 for (UInt_t irow=0;irow<nrows;++irow){
334 UInt_t npads=noiseROC->GetNPads(irow);
335 for (UInt_t ipad=0;ipad<npads;++ipad){
336 //don't use masked channels;
337 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
338 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
339 //check if noise==0
340 if (noiseVal==0) {
341 ++nonMaskedZero;
342 continue;
343 }
344 //check for nan
345 if ( !(noiseVal<10000000) ){
346 printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
347 continue;
348 }
349 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
350 Int_t masksen=1; // sensitive pards are not masked (0)
351 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
352 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
353 //IROCs
354 if (irow>19&&irow<46){
355 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
356 }
357 Int_t type=1;
358 vNoiseMean[type]+=noiseVal;
359 vNoiseRMS[type]+=noiseVal*noiseVal;
360 ++c[type];
361 if (!masksen){
362 vNoiseMeanSenRegions[type]+=noiseVal;
363 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
364 ++cs[type];
365 }
366 } else {
367 //OROCs
368 //define sensive regions
369 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
370 if ( irow>75 ){
371 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
372 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
373 }
374 if ((Int_t)irow<par.GetNRowUp1()){
375 //OROC1
376 Int_t type=2;
377 vNoiseMean[type]+=noiseVal;
378 vNoiseRMS[type]+=noiseVal*noiseVal;
379 ++c[type];
380 if (!masksen){
381 vNoiseMeanSenRegions[type]+=noiseVal;
382 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
383 ++cs[type];
384 }
385 }else{
386 //OROC2
387 Int_t type=3;
388 vNoiseMean[type]+=noiseVal;
389 vNoiseRMS[type]+=noiseVal*noiseVal;
390 ++c[type];
391 if (!masksen){
392 vNoiseMeanSenRegions[type]+=noiseVal;
393 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
394 ++cs[type];
395 }
396 }
397 }
398 //whole tpc
399 Int_t type=0;
400 vNoiseMean[type]+=noiseVal;
401 vNoiseRMS[type]+=noiseVal*noiseVal;
402 ++c[type];
403 if (!masksen){
404 vNoiseMeanSenRegions[type]+=noiseVal;
405 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
406 ++cs[type];
407 }
408 }//end loop pads
409 }//end loop rows
410 }//end loop sectors (rocs)
411
412 //calculate mean and RMS
413 const Double_t verySmall=0.0000000001;
414 for (UInt_t i=0;i<infoSize;++i){
415 Double_t mean=0;
416 Double_t rms=0;
417 Double_t meanSen=0;
418 Double_t rmsSen=0;
419
420 if (c[i]>verySmall){
421// printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
422 mean=vNoiseMean[i]/c[i];
423 rms=vNoiseRMS[i];
424 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
425 }
426 vNoiseMean[i]=mean;
427 vNoiseRMS[i]=rms;
428
429 if (cs[i]>verySmall){
430 meanSen=vNoiseMeanSenRegions[i]/cs[i];
431 rmsSen=vNoiseRMSSenRegions[i];
432 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
433 }
434 vNoiseMeanSenRegions[i]=meanSen;
435 vNoiseRMSSenRegions[i]=rmsSen;
436 }
437}
438
439//_____________________________________________________________________________________
440void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
441{
442 //
443 // Process the Pulser information
444 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
445 //
446
447 const UInt_t infoSize=4;
448 //reset counters to error number
449 vMeanTime.ResizeTo(infoSize);
450 vMeanTime.Zero();
451 //counter
452 TVectorD c(infoSize);
453 //retrieve pulser and ALTRO data
454 if (!fPulserTmean) return;
455 //
456 //get Outliers
457 AliTPCCalROC *rocOut=0x0;
458 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
459 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
460 if (!tmeanROC) continue;
461 rocOut=fPulserOutlier->GetCalROC(isec);
462 UInt_t nchannels=tmeanROC->GetNchannels();
463 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
464 if (rocOut && rocOut->GetValue(ichannel)) continue;
465 Float_t val=tmeanROC->GetValue(ichannel);
466 Int_t type=isec/18;
467 vMeanTime[type]+=val;
468 ++c[type];
469 }
470 }
471 //calculate mean
472 for (UInt_t itype=0; itype<infoSize; ++itype){
473 if (c[itype]>0) vMeanTime[itype]/=c[itype];
474 else vMeanTime[itype]=0;
475 }
476}
477//_____________________________________________________________________________________
478void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
479{
480 //
481 // Get Values from ALTRO configuration data
482 //
483 nMasked=-1;
484 if (!fALTROMasked) return;
485 nMasked=0;
486 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
487 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
488 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
489 if (rocMasked->GetValue(ichannel)) ++nMasked;
490 }
491 }
492}
493//_____________________________________________________________________________________
494void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
495{
496 //
497 // Proces Goofie values, return statistical information of the currently set goofieArray
498 // The meaning of the entries are given below
499 /*
500 1 TPC_ANODE_I_A00_STAT
501 2 TPC_DVM_CO2
502 3 TPC_DVM_DriftVelocity
503 4 TPC_DVM_FCageHV
504 5 TPC_DVM_GainFar
505 6 TPC_DVM_GainNear
506 7 TPC_DVM_N2
507 8 TPC_DVM_NumberOfSparks
508 9 TPC_DVM_PeakAreaFar
509 10 TPC_DVM_PeakAreaNear
510 11 TPC_DVM_PeakPosFar
511 12 TPC_DVM_PeakPosNear
512 13 TPC_DVM_PickupHV
513 14 TPC_DVM_Pressure
514 15 TPC_DVM_T1_Over_P
515 16 TPC_DVM_T2_Over_P
516 17 TPC_DVM_T_Over_P
517 18 TPC_DVM_TemperatureS1
518 */
519 if (!fGoofieArray){
520 Int_t nsensors=19;
521 vecEntries.ResizeTo(nsensors);
522 vecMedian.ResizeTo(nsensors);
523 vecMean.ResizeTo(nsensors);
524 vecRMS.ResizeTo(nsensors);
525 vecEntries.Zero();
526 vecMedian.Zero();
527 vecMean.Zero();
528 vecRMS.Zero();
529 return;
530 }
531 Double_t kEpsilon=0.0000000001;
532 Double_t kBig=100000000000.;
533 Int_t nsensors = fGoofieArray->NumSensors();
534 vecEntries.ResizeTo(nsensors);
535 vecMedian.ResizeTo(nsensors);
536 vecMean.ResizeTo(nsensors);
537 vecRMS.ResizeTo(nsensors);
538 TVectorF values;
539 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
540 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
541 if (gsensor && gsensor->GetGraph()){
542 Int_t npoints = gsensor->GetGraph()->GetN();
543 // filter zeroes
544 values.ResizeTo(npoints);
545 Int_t nused =0;
546 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
547 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
548 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
549 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
550 nused++;
551 }
552 }
553 //
554 vecEntries[isensor]= nused;
555 if (nused>1){
556 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
557 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
558 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
559 }
560 }
561 }
562}
7390f655 563//_____________________________________________________________________________________
564void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
565{
566 //
567 // check the variations of the pedestal data to the reference pedestal data
568 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
569 //
570 const Int_t npar=4;
571 TVectorF vThres(npar); //thresholds
572 Int_t nActive=0; //number of active channels
573
574 //reset and set thresholds
575 pedestalDeviations.ResizeTo(npar);
576 for (Int_t i=0;i<npar;++i){
577 pedestalDeviations.GetMatrixArray()[i]=0;
578 vThres.GetMatrixArray()[i]=(i+1)*.5;
579 }
580 //check all needed data is available
581 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
582 //loop over all channels
583 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
584 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
585 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
586 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
587 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
588 UInt_t nrows=mROC->GetNrows();
589 for (UInt_t irow=0;irow<nrows;++irow){
590 UInt_t npads=mROC->GetNPads(irow);
591 for (UInt_t ipad=0;ipad<npads;++ipad){
592 //don't use masked channels;
593 if (mROC ->GetValue(irow,ipad)) continue;
594 if (mRefROC->GetValue(irow,ipad)) continue;
595 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
596 for (Int_t i=0;i<npar;++i){
597 if (deviation>vThres[i])
598 ++pedestalDeviations.GetMatrixArray()[i];
599 }
600 ++nActive;
601 }//end ipad
602 }//ind irow
603 }//end isec
604 if (nActive>0){
605 for (Int_t i=0;i<npar;++i){
606 pedestalDeviations.GetMatrixArray()[i]/=nActive;
607 }
608 }
609}
610//_____________________________________________________________________________________
611void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
612{
613 //
614 // check the variations of the noise data to the reference noise data
615 // thresholds are 5, 10, 15 and 20 percent respectively.
616 //
617 const Int_t npar=4;
618 TVectorF vThres(npar); //thresholds
619 Int_t nActive=0; //number of active channels
620
621 //reset and set thresholds
622 noiseDeviations.ResizeTo(npar);
623 for (Int_t i=0;i<npar;++i){
624 noiseDeviations.GetMatrixArray()[i]=0;
625 vThres.GetMatrixArray()[i]=(i+1)*.05;
626 }
627 //check all needed data is available
628 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
629 //loop over all channels
630 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
631 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
632 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
633 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
634 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
635 UInt_t nrows=mROC->GetNrows();
636 for (UInt_t irow=0;irow<nrows;++irow){
637 UInt_t npads=mROC->GetNPads(irow);
638 for (UInt_t ipad=0;ipad<npads;++ipad){
639 //don't use masked channels;
640 if (mROC ->GetValue(irow,ipad)) continue;
641 if (mRefROC->GetValue(irow,ipad)) continue;
642 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
643 for (Int_t i=0;i<npar;++i){
644 if (deviation>vThres[i])
645 ++noiseDeviations.GetMatrixArray()[i];
646 }
647 ++nActive;
648 }//end ipad
649 }//ind irow
650 }//end isec
651 if (nActive>0){
652 for (Int_t i=0;i<npar;++i){
653 noiseDeviations.GetMatrixArray()[i]/=nActive;
654 }
655 }
656}
657//_____________________________________________________________________________________
658void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
659 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
660{
661 //
662 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
663 // thresholds are .5, 1, 5 and 10 percent respectively.
664 //
665 //
666 const Int_t npar=4;
667 TVectorF vThres(npar); //thresholds
668 Int_t nActive=0; //number of active channels
669
670 //reset and set thresholds
671 pulserQdeviations.ResizeTo(npar);
672 for (Int_t i=0;i<npar;++i){
673 pulserQdeviations.GetMatrixArray()[i]=0;
674 }
675 npadsOutOneTB=0;
676 npadsOffAdd=0;
677 varQMean=0;
678 vThres.GetMatrixArray()[0]=.005;
679 vThres.GetMatrixArray()[1]=.01;
680 vThres.GetMatrixArray()[2]=.05;
681 vThres.GetMatrixArray()[3]=.1;
682 //check all needed data is available
683 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
684 //
685 UpdateRefPulserOutlierMap();
686 //loop over all channels
732e90a8 687 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 688 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
689 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
690 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
732e90a8 691// AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
7390f655 692 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
693 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
694 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
695 Float_t pt_mean=ptROC->GetMean(oROC);
696 UInt_t nrows=mROC->GetNrows();
697 for (UInt_t irow=0;irow<nrows;++irow){
698 UInt_t npads=mROC->GetNPads(irow);
699 for (UInt_t ipad=0;ipad<npads;++ipad){
700 //don't use masked channels;
701 if (mROC ->GetValue(irow,ipad)) continue;
702 if (mRefROC->GetValue(irow,ipad)) continue;
703 //don't user edge pads
704 if (ipad==0||ipad==npads-1) continue;
705 //data
706 Float_t pq=pqROC->GetValue(irow,ipad);
707 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
708 Float_t pt=ptROC->GetValue(irow,ipad);
709// Float_t ptRef=ptRefROC->GetValue(irow,ipad);
710 //comparisons q
711 Float_t deviation=TMath::Abs(pq/pqRef-1);
712 for (Int_t i=0;i<npar;++i){
713 if (deviation>vThres[i])
714 ++pulserQdeviations.GetMatrixArray()[i];
715 }
716 if (pqRef>11&&pq<11) ++npadsOffAdd;
717 varQMean+=pq-pqRef;
718 //comparisons t
719 if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB;
720 ++nActive;
721 }//end ipad
722 }//ind irow
723 }//end isec
724 if (nActive>0){
725 for (Int_t i=0;i<npar;++i){
726 pulserQdeviations.GetMatrixArray()[i]/=nActive;
727 varQMean/=nActive;
728 }
729 }
730}
892226be 731//_____________________________________________________________________________________
732void AliTPCcalibDButil::UpdatePulserOutlierMap()
7390f655 733{
734 //
735 //
736 //
737 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
738}
739//_____________________________________________________________________________________
740void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
741{
742 //
743 //
744 //
745 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
746}
747//_____________________________________________________________________________________
748void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
892226be 749{
750 //
751 // Create a map that contains outliers from the Pulser calibration data.
752 // The outliers include masked channels, edge pads and pads with
753 // too large timing and charge variations.
7390f655 754 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
892226be 755 // those do not contain masked and edge pads
756 //
7390f655 757 if (!pulT||!pulQ) {
892226be 758 //reset map
7390f655 759 pulOut->Multiply(0.);
892226be 760 fNpulserOutliers=-1;
761 return;
762 }
763 AliTPCCalROC *rocMasked=0x0;
764 fNpulserOutliers=0;
765
766 //Create Outlier Map
767 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 768 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
769 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
770 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
892226be 771 if (!tmeanROC||!qmeanROC) {
772 //reset outliers in this ROC
773 outROC->Multiply(0.);
774 continue;
775 }
776 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
777// Double_t dummy=0;
778// Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
779// Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
780 UInt_t nrows=tmeanROC->GetNrows();
781 for (UInt_t irow=0;irow<nrows;++irow){
782 UInt_t npads=tmeanROC->GetNPads(irow);
783 for (UInt_t ipad=0;ipad<npads;++ipad){
784 Int_t outlier=0,masked=0;
785 Float_t q=qmeanROC->GetValue(irow,ipad);
786 Float_t t=tmeanROC->GetValue(irow,ipad);
787 //masked channels are outliers
788 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
789 //edge pads are outliers
790 if (ipad==0||ipad==npads-1) masked=1;
791 //channels with too large charge or timing deviation from the meadian are outliers
792// if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
793 if (q<fPulQminLimit && !masked) outlier=1;
794 //check for nan
795 if ( !(q<10000000) || !(t<10000000)) outlier=1;
796 outROC->SetValue(irow,ipad,outlier+masked);
797 fNpulserOutliers+=outlier;
798 }
799 }
800 }
801}
802//_____________________________________________________________________________________
2cb269df 803AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
892226be 804{
805 //
806 // Create pad time0 object from pulser and/or CE data, depending on the selected model
807 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
808 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
732e90a8 809 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
892226be 810 //
2cb269df 811 // In case model 2 is invoked - gy arival time gradient is also returned
812 //
813 gyA=0;
814 gyC=0;
892226be 815 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
816 // decide between different models
817 if (model==0||model==1){
818 TVectorD vMean;
819 if (model==1) ProcessPulser(vMean);
820 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
821 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
822 if (!rocPulTmean) continue;
823 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
824 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
825 Float_t mean=rocPulTmean->GetMean(rocOut);
826 //treat case where a whole partition is masked
827 if (mean==0) mean=rocPulTmean->GetMean();
828 if (model==1) {
829 Int_t type=isec/18;
830 mean=vMean[type];
831 }
832 UInt_t nrows=rocTime0->GetNrows();
833 for (UInt_t irow=0;irow<nrows;++irow){
834 UInt_t npads=rocTime0->GetNPads(irow);
835 for (UInt_t ipad=0;ipad<npads;++ipad){
836 Float_t time=rocPulTmean->GetValue(irow,ipad);
837 //in case of an outlier pad use the mean of the altro values.
838 //This should be the most precise guess in that case.
839 if (rocOut->GetValue(irow,ipad)) {
840 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
841 if (time==0) time=mean;
842 }
843 Float_t val=time-mean;
844 rocTime0->SetValue(irow,ipad,val);
845 }
846 }
847 }
2cb269df 848 } else if (model==2){
849 Double_t pgya,pgyc,pchi2a,pchi2c;
850 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
851 fCETmean->Add(padPulser,-1.);
732e90a8 852 TVectorD vA,vC;
853 AliTPCCalPad outCE("outCE","outCE");
854 Int_t nOut;
2cb269df 855 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
856 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
732e90a8 857// AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
2cb269df 858 if (!padFit) { delete padPulser; return 0;}
859 gyA=vA[2];
860 gyC=vC[2];
861 fCETmean->Add(padPulser,1.);
732e90a8 862 padTime0->Add(fCETmean);
2cb269df 863 padTime0->Add(padFit,-1);
864 delete padPulser;
732e90a8 865 TVectorD vFitROC;
866 TMatrixD mFitROC;
867 Float_t chi2;
868 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
869 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
870 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
871 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
872 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
873 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
874 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
875 Float_t mean=rocPulTmean->GetMean(rocOutPul);
876 if (mean==0) mean=rocPulTmean->GetMean();
877 UInt_t nrows=rocTime0->GetNrows();
878 for (UInt_t irow=0;irow<nrows;++irow){
879 UInt_t npads=rocTime0->GetNPads(irow);
880 for (UInt_t ipad=0;ipad<npads;++ipad){
881 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
882 if (rocOutCE->GetValue(irow,ipad)){
883 Float_t valOut=rocCEfit->GetValue(irow,ipad);
884 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
885 rocTime0->SetValue(irow,ipad,valOut);
886 }
887 }
888 }
889 delete rocCEfit;
890 }
891 delete padFit;
892226be 892 }
2cb269df 893 Double_t median = padTime0->GetMedian();
894 padTime0->Add(-median); // normalize to median
892226be 895 return padTime0;
896}
897//_____________________________________________________________________________________
898Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
899{
900 if (roc==0) return 0.;
901 const Int_t sector=roc->GetSector();
902 AliTPCROC *tpcRoc=AliTPCROC::Instance();
903 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
904 Float_t mean=0;
905 Int_t n=0;
906
907 //loop over a small range around the requested pad (+-10 rows/pads)
908 for (Int_t irow=row-10;irow<row+10;++irow){
909 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
910 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
911 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
912 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
913 if (altroRoc!=altroCurr) continue;
914 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
915 Float_t val=roc->GetValue(irow,ipad);
916 mean+=val;
917 ++n;
918 }
919 }
920 if (n>0) mean/=n;
921 return mean;
922}
7390f655 923//_____________________________________________________________________________________
924void AliTPCcalibDButil::SetRefFile(const char* filename)
925{
926 //
927 // load cal pad objects form the reference file
928 //
929 TDirectory *currDir=gDirectory;
930 TFile f(filename);
931 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
932 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
933 //pulser data
934 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
935 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
936 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
937 //CE data
938 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
939 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
940 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
941 //Altro data
942// fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
943// fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
944// fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
945// fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
946 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
947 f.Close();
948 currDir->cd();
949}
950
2cb269df 951
952
953
954AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad *ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
955 //
956 // Author: marian.ivanov@cern.ch
957 //
958 // Create outlier map for CE study
959 // Parameters:
960 // Return value - outlyer map
961 // noutlyersCE - number of outlyers
962 // minSignal - minimal total Q signal
963 // cutRMSMin - minimal width of the signal in respect to the median
964 // cutRMSMax - maximal width of the signal in respect to the median
965 // cutMaxDistT - maximal deviation from time median per chamber
966 //
967 // Outlyers criteria:
968 // 0. Exclude masked pads
969 // 1. Exclude first two rows in IROC and last two rows in OROC
970 // 2. Exclude edge pads
971 // 3. Exclude channels with too large variations
972 // 4. Exclude pads with too small signal
973 // 5. Exclude signal with outlyers RMS
974 // 6. Exclude channels to far from the chamber median
975 noutliersCE=0;
976 //create outlier map
977 AliTPCCalPad *out=ceOut;
978 if (!out) out= new AliTPCCalPad("outCE","outCE");
979 AliTPCCalROC *rocMasked=0x0;
980 if (!fCETmean) return 0;
981 if (!fCETrms) return 0;
982 if (!fCEQmean) return 0;
983 //
984 //loop over all channels
985 //
986 Double_t rmsMedian = fCETrms->GetMedian();
987 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
988 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
989 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
990 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
991 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
992 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
993 Double_t trocMedian = rocData->GetMedian();
994 //
995 if (!rocData) {
996 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
997 rocOut->Add(1.);
998 continue;
999 }
1000 //
1001 //select outliers
1002 UInt_t nrows=rocData->GetNrows();
1003 for (UInt_t irow=0;irow<nrows;++irow){
1004 UInt_t npads=rocData->GetNPads(irow);
1005 for (UInt_t ipad=0;ipad<npads;++ipad){
1006 rocOut->SetValue(irow,ipad,0);
1007 Float_t valTmean=rocData->GetValue(irow,ipad);
1008 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1009 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1010 //0. exclude masked pads
1011 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1012 rocOut->SetValue(irow,ipad,1);
1013 continue;
1014 }
1015 //1. exclude first two rows in IROC and last two rows in OROC
1016 if (iroc<36){
1017 if (irow<2) rocOut->SetValue(irow,ipad,1);
1018 } else {
1019 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1020 }
1021 //2. exclude edge pads
1022 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1023 //exclude values that are exactly 0
1024 if (valTmean==0) {
1025 rocOut->SetValue(irow,ipad,1);
1026 ++noutliersCE;
1027 }
1028 //3. exclude channels with too large variations
1029 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1030 rocOut->SetValue(irow,ipad,1);
1031 ++noutliersCE;
1032 }
1033 //
1034 //4. exclude channels with too small signal
1035 if (valQmean<minSignal) {
1036 rocOut->SetValue(irow,ipad,1);
1037 ++noutliersCE;
1038 }
1039 //
1040 //5. exclude channels with too small rms
1041 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1042 rocOut->SetValue(irow,ipad,1);
1043 ++noutliersCE;
1044 }
1045 //
1046 //6. exclude channels to far from the chamber median
1047 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1048 rocOut->SetValue(irow,ipad,1);
1049 ++noutliersCE;
1050 }
1051 }
1052 }
1053 }
1054 //
1055 return out;
1056}
1057
1058
1059AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad *pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1060 //
1061 // Author: marian.ivanov@cern.ch
1062 //
1063 // Create outlier map for Pulser
1064 // Parameters:
1065 // Return value - outlyer map
1066 // noutlyersPulser - number of outlyers
1067 // cutTime - absolute cut - distance to the median of chamber
1068 // cutnRMSQ - nsigma cut from median q distribution per chamber
1069 // cutnRMSrms - nsigma cut from median rms distribution
1070 // Outlyers criteria:
1071 // 0. Exclude masked pads
1072 // 1. Exclude time outlyers (default 3 time bins)
1073 // 2. Exclude q outlyers (default 5 sigma)
1074 // 3. Exclude rms outlyers (default 5 sigma)
1075 noutliersPulser=0;
1076 AliTPCCalPad *out=pulserOut;
1077 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1078 AliTPCCalROC *rocMasked=0x0;
1079 if (!fPulserTmean) return 0;
1080 if (!fPulserTrms) return 0;
1081 if (!fPulserQmean) return 0;
1082 //
1083 //loop over all channels
1084 //
1085 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1086 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1087 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1088 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1089 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1090 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1091 //
1092 Double_t rocMedianT = rocData->GetMedian();
1093 Double_t rocMedianQ = rocPulserQ->GetMedian();
1094 Double_t rocRMSQ = rocPulserQ->GetRMS();
1095 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1096 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1097 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1098 rocOut->SetValue(ichannel,0);
1099 Float_t valTmean=rocData->GetValue(ichannel);
1100 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1101 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1102 Int_t isOut=0;
1103 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1104 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1105 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1106 rocOut->SetValue(ichannel,isOut);
1107 if (isOut) noutliersPulser++;
1108 }
1109 }
1110 return out;
1111}
1112
1113
1114AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1115 //
1116 // Author : Marian Ivanov
1117 // Create pad time0 correction map using information from the CE and from pulser
1118 //
1119 //
1120 // Return PadTime0 to be used for time0 relative alignment
1121 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1122 // using $ALICE_ROOT/TPC/script/gui application
1123 //
1124 // fitResultsA - fitParameters A side
1125 // fitResultsC - fitParameters C side
1126 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1127 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1128 //
1129 //
1130 // Algorithm:
1131 // 1. Find outlier map for CE
1132 // 2. Find outlier map for Pulser
1133 // 3. Replace outlier by median at given sector (median without outliers)
1134 // 4. Substract from the CE data pulser
1135 // 5. Fit the CE with formula
1136 // 5.1) (IROC-OROC) offset
1137 // 5.2) gx
1138 // 5.3) gy
1139 // 5.4) (lx-xmid)
1140 // 5.5) (IROC-OROC)*(lx-xmid)
1141 // 5.6) (ly/lx)^2
1142 // 6. Substract gy fit dependence from the CE data
1143 // 7. Add pulser back to CE data
1144 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1145 // 9. return CE data
1146 //
1147 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1148 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1149
1150 // fit formula
1151 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)";
1152 // output for fit formula
1153 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)";
1154 // gy part of formula
1155 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)";
1156 //
1157 //
1158 if (!fCETmean) return 0;
1159 Double_t pgya,pgyc,pchi2a,pchi2c;
1160 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1161 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1162
1163 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1164 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1165 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1166 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1167 padPulser->SetName("padPulser");
1168 padPulserOut->SetName("padPulserOut");
1169 padCE->SetName("padCE");
1170 padCEIn->SetName("padCEIn");
1171 padCEOut->SetName("padCEOut");
1172 padOut->SetName("padOut");
1173
1174 //
1175 // make combined outlyers map
1176 // and replace outlyers in maps with median for chamber
1177 //
1178 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1179 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1180 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1181 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1182 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1183 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1184 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1185 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1186 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1187 if (rocPulserOut->GetValue(ichannel)>0) {
1188 rocPulser->SetValue(ichannel,pulserMedian);
1189 rocOut->SetValue(ichannel,1);
1190 }
1191 if (rocCEOut->GetValue(ichannel)>0) {
1192 rocCE->SetValue(ichannel,ceMedian);
1193 rocOut->SetValue(ichannel,1);
1194 }
1195 }
1196 }
1197 //
1198 // remove pulser time 0
1199 //
1200 padCE->Add(padPulser,-1);
1201 //
1202 // Make fits
1203 //
1204 TMatrixD dummy;
1205 Float_t chi2Af,chi2Cf;
1206 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1207 chi2A=chi2Af;
1208 chi2C=chi2Cf;
1209 //
1210 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1211 padCEFitGY->SetName("padCEFitGy");
1212 //
1213 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1214 padCEFit->SetName("padCEFit");
1215 //
1216 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1217 padCEDiff->SetName("padCEDiff");
1218 padCEDiff->Add(padCEFit,-1.);
1219 //
1220 //
1221 padCE->Add(padCEFitGY,-1.);
1222
1223 padCE->Add(padPulser,1.);
1224 Double_t padmedian = padCE->GetMedian();
1225 padCE->Add(-padmedian); // normalize to median
1226 //
1227 // Replace outliers by fit value - median of diff per given chamber -GY fit
1228 //
1229 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1230 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1231 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1232 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1233 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1234 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1235 //
1236 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1237 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1238 if (rocOut->GetValue(ichannel)==0) continue;
1239 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1240 rocCE->SetValue(ichannel,value);
1241 }
1242 }
1243 //
1244 //
1245 if (dumpfile){
1246 //dump to the file - result can be visualized
1247 AliTPCPreprocessorOnline preprocesor;
1248 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1249 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1250 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1251 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1252 //
1253 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1254 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1255 //
1256 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1257 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1258 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1259 preprocesor.DumpToFile(dumpfile);
1260 }
1261 delete padPulser;
1262 delete padPulserOut;
1263 delete padCEIn;
1264 delete padCEOut;
1265 delete padOut;
1266 delete padCEDiff;
1267 delete padCEFitGY;
1268 return padCE;
1269}
1270
817766d5 1271
1272
1273
1274
1275Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1276 //
1277 // find the closest point to xref in x direction
1278 // return dx and value
1279 Int_t index=0;
1280 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1281 if (index<0) index=0;
1282 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1283 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1284 dx = xref-graph->GetX()[index];
1285 y = graph->GetY()[index];
1286 return index;
1287}
1288
1289
1290Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1291 //
1292 // Get the correction of the trigger offset
1293 // combining information from the laser track calibration
1294 // and from cosmic calibration
1295 //
1296 // run - run number
1297 // timeStamp - tim stamp in seconds
1298 // deltaT - integration period to calculate offset
1299 // deltaTLaser -max validity of laser data
1300 // valType - 0 - median, 1- mean
1301 //
1302 // Integration vaues are just recomendation - if not possible to get points
1303 // automatically increase the validity by factor 2
1304 // (recursive algorithm until one month of data taking)
1305 //
1306 //
1307 const Float_t kLaserCut=0.0005;
1308 const Int_t kMaxPeriod=3600*24*30*3; // 3 month max
1309 const Int_t kMinPoints=20;
1310 //
1311 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1312 if (!array) {
1313 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1314 }
1315 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1316 if (!array) return 0;
1317 //
1318 TGraphErrors *laserA[3]={0,0,0};
1319 TGraphErrors *laserC[3]={0,0,0};
1320 TGraphErrors *cosmicAll=0;
1321 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1322 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1323 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1324 //
1325 //
1326 if (!cosmicAll) return 0;
1327 Int_t nmeasC=cosmicAll->GetN();
1328 Float_t *tdelta = new Float_t[nmeasC];
1329 Int_t nused=0;
1330 for (Int_t i=0;i<nmeasC;i++){
1331 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1332 Float_t ccosmic=cosmicAll->GetY()[i];
1333 Double_t yA=0,yC=0,dA=0,dC=0;
1334 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1335 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1336 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1337 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1338 //
1339 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1340 Float_t claser=0;
1341 if (TMath::Abs(yA-yC)<kLaserCut) {
1342 claser=(yA-yC)*0.5;
1343 }else{
1344 if (i%2==0) claser=yA;
1345 if (i%2==1) claser=yC;
1346 }
1347 tdelta[nused]=ccosmic-claser;
1348 nused++;
1349 }
1350 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1351 Double_t median = TMath::Median(nused,tdelta);
1352 Double_t mean = TMath::Mean(nused,tdelta);
1353 delete tdelta;
1354 return (valType==0) ? median:mean;
1355}
1356
1357Double_t AliTPCcalibDButil::GetVDriftTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1358 //
1359 // Get the correction of the drift velocity
1360 // combining information from the laser track calibration
1361 // and from cosmic calibration
1362 //
1363 // run - run number
1364 // timeStamp - tim stamp in seconds
1365 // deltaT - integration period to calculate time0 offset
1366 // deltaTLaser -max validity of laser data
1367 // valType - 0 - median, 1- mean
1368 //
1369 // Integration vaues are just recomendation - if not possible to get points
1370 // automatically increase the validity by factor 2
1371 // (recursive algorithm until one month of data taking)
1372 //
1373 //
1374 //
1375 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1376 if (!array) {
1377 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1378 }
1379 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1380 if (!array) return 0;
1381 TGraphErrors *cosmicAll=0;
1382 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1383 if (!cosmicAll) return 0;
1384 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1385 Double_t vcosmic= cosmicAll->Eval(timeStamp);
1386 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=timeStamp>cosmicAll->GetY()[cosmicAll->GetN()-1];
1387 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1388 return vcosmic+t0;
1389
1390 /*
1391 Example usage:
1392
1393 Int_t run=89000
1394 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1395 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1396 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1397 //
1398 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1399 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1400 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1401 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1402
1403 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1404 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1405
1406 */
1407
1408}
1409