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