ATO-78 a.) adding aliases used for gain correction fitting \n b)Kr and QA prefix...
[u/mrichter/AliRoot.git] / TPC / Base / 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>
949d8707 33#include <TMap.h>
34#include <TGraphErrors.h>
a23ba1c3 35#include <AliCDBStorage.h>
892226be 36#include <AliDCSSensorArray.h>
a23ba1c3 37#include <AliTPCSensorTempArray.h>
892226be 38#include <AliDCSSensor.h>
949d8707 39#include <AliLog.h>
40#include <AliCDBEntry.h>
41#include <AliCDBManager.h>
42#include <AliCDBId.h>
df5d3222 43#include <AliSplineFit.h>
892226be 44#include "AliTPCcalibDB.h"
45#include "AliTPCCalPad.h"
46#include "AliTPCCalROC.h"
47#include "AliTPCROC.h"
48#include "AliTPCmapper.h"
49#include "AliTPCParam.h"
6e7d7dc4 50#include "AliTPCCalibRaw.h"
949d8707 51#include "AliTPCPreprocessorOnline.h"
52#include "AliTPCdataQA.h"
a23ba1c3 53#include "AliLog.h"
892226be 54#include "AliTPCcalibDButil.h"
a23ba1c3 55#include "AliTPCCalibVdrift.h"
1fabc823 56#include "AliMathBase.h"
a980538f 57#include "AliRelAlignerKalman.h"
c4e4498c 58#include "TTree.h"
59#include "TROOT.h"
c4f17a0d 60#include "TKey.h"
892226be 61
8166d5d7 62const Float_t kAlmost0=1.e-30;
63
892226be 64ClassImp(AliTPCcalibDButil)
65AliTPCcalibDButil::AliTPCcalibDButil() :
66 TObject(),
1e722a63 67 fCalibDB(0),
892226be 68 fPadNoise(0x0),
69 fPedestals(0x0),
70 fPulserTmean(0x0),
71 fPulserTrms(0x0),
72 fPulserQmean(0x0),
73 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
74 fCETmean(0x0),
75 fCETrms(0x0),
76 fCEQmean(0x0),
77 fALTROMasked(0x0),
6e7d7dc4 78 fCalibRaw(0x0),
949d8707 79 fDataQA(0x0),
80 fRefMap(0x0),
81 fCurrentRefMap(0x0),
82 fRefValidity(""),
7390f655 83 fRefPadNoise(0x0),
84 fRefPedestals(0x0),
949d8707 85 fRefPedestalMasked(0x0),
7390f655 86 fRefPulserTmean(0x0),
87 fRefPulserTrms(0x0),
88 fRefPulserQmean(0x0),
89 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
949d8707 90 fRefPulserMasked(0x0),
7390f655 91 fRefCETmean(0x0),
92 fRefCETrms(0x0),
93 fRefCEQmean(0x0),
949d8707 94 fRefCEMasked(0x0),
95 fRefALTROFPED(0x0),
96 fRefALTROZsThr(0x0),
97 fRefALTROAcqStart(0x0),
98 fRefALTROAcqStop(0x0),
7390f655 99 fRefALTROMasked(0x0),
100 fRefCalibRaw(0x0),
949d8707 101 fRefDataQA(0x0),
892226be 102 fGoofieArray(0x0),
103 fMapper(new AliTPCmapper(0x0)),
104 fNpulserOutliers(-1),
7390f655 105 fIrocTimeOffset(0),
892226be 106 fCETmaxLimitAbs(1.5),
107 fPulTmaxLimitAbs(1.5),
108 fPulQmaxLimitAbs(5),
a23ba1c3 109 fPulQminLimit(11),
110 fRuns(0), // run list with OCDB info
111 fRunsStart(0), // start time for given run
112 fRunsStop(0) // stop time for given run
892226be 113{
114 //
115 // Default ctor
116 //
117}
118//_____________________________________________________________________________________
119AliTPCcalibDButil::~AliTPCcalibDButil()
120{
121 //
122 // dtor
123 //
124 delete fPulserOutlier;
7390f655 125 delete fRefPulserOutlier;
892226be 126 delete fMapper;
7390f655 127 if (fRefPadNoise) delete fRefPadNoise;
128 if (fRefPedestals) delete fRefPedestals;
949d8707 129 if (fRefPedestalMasked) delete fRefPedestalMasked;
7390f655 130 if (fRefPulserTmean) delete fRefPulserTmean;
131 if (fRefPulserTrms) delete fRefPulserTrms;
132 if (fRefPulserQmean) delete fRefPulserQmean;
949d8707 133 if (fRefPulserMasked) delete fRefPulserMasked;
7390f655 134 if (fRefCETmean) delete fRefCETmean;
135 if (fRefCETrms) delete fRefCETrms;
136 if (fRefCEQmean) delete fRefCEQmean;
949d8707 137 if (fRefCEMasked) delete fRefCEMasked;
138 if (fRefALTROFPED) delete fRefALTROFPED;
139 if (fRefALTROZsThr) delete fRefALTROZsThr;
140 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
141 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
7390f655 142 if (fRefALTROMasked) delete fRefALTROMasked;
143 if (fRefCalibRaw) delete fRefCalibRaw;
949d8707 144 if (fCurrentRefMap) delete fCurrentRefMap;
892226be 145}
146//_____________________________________________________________________________________
147void AliTPCcalibDButil::UpdateFromCalibDB()
148{
149 //
150 // Update pointers from calibDB
151 //
1e722a63 152 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
56ce896d 153 fCalibDB->UpdateNonRec(); // load all infromation now
892226be 154 fPadNoise=fCalibDB->GetPadNoise();
155 fPedestals=fCalibDB->GetPedestals();
156 fPulserTmean=fCalibDB->GetPulserTmean();
157 fPulserTrms=fCalibDB->GetPulserTrms();
158 fPulserQmean=fCalibDB->GetPulserQmean();
159 fCETmean=fCalibDB->GetCETmean();
160 fCETrms=fCalibDB->GetCETrms();
161 fCEQmean=fCalibDB->GetCEQmean();
162 fALTROMasked=fCalibDB->GetALTROMasked();
163 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
6e7d7dc4 164 fCalibRaw=fCalibDB->GetCalibRaw();
949d8707 165 fDataQA=fCalibDB->GetDataQA();
892226be 166 UpdatePulserOutlierMap();
22c558aa 167// SetReferenceRun();
168// UpdateRefDataFromOCDB();
892226be 169}
170//_____________________________________________________________________________________
7390f655 171void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
8166d5d7 172 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
892226be 173{
174 //
175 // Process the CE data for this run
176 // the return TVectorD arrays contian the results of the fit
177 // noutliersCE contains the number of pads marked as outliers,
178 // not including masked and edge pads
179 //
180
181 //retrieve CE and ALTRO data
182 if (!fCETmean){
183 TString fitString(fitFormula);
184 fitString.ReplaceAll("++","#");
185 Int_t ndim=fitString.CountChar('#')+2;
186 fitResultsA.ResizeTo(ndim);
187 fitResultsC.ResizeTo(ndim);
188 fitResultsA.Zero();
189 fitResultsC.Zero();
190 noutliersCE=-1;
191 return;
192 }
193 noutliersCE=0;
194 //create outlier map
7390f655 195 AliTPCCalPad *out=0;
196 if (outCE) out=outCE;
197 else out=new AliTPCCalPad("outCE","outCE");
892226be 198 AliTPCCalROC *rocMasked=0x0;
199 //loop over all channels
200 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
201 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
202 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
7390f655 203 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
6e7d7dc4 204 if (!rocData) {
205 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
7390f655 206 rocOut->Add(1.);
6e7d7dc4 207 continue;
208 }
892226be 209 //add time offset to IROCs
210 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
211 rocData->Add(fIrocTimeOffset);
212 //select outliers
213 UInt_t nrows=rocData->GetNrows();
214 for (UInt_t irow=0;irow<nrows;++irow){
215 UInt_t npads=rocData->GetNPads(irow);
216 for (UInt_t ipad=0;ipad<npads;++ipad){
7390f655 217 rocOut->SetValue(irow,ipad,0);
892226be 218 //exclude masked pads
219 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
7390f655 220 rocOut->SetValue(irow,ipad,1);
892226be 221 continue;
222 }
732e90a8 223 //exclude first two rows in IROC and last two rows in OROC
224 if (iroc<36){
225 if (irow<2) rocOut->SetValue(irow,ipad,1);
226 } else {
227 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
228 }
892226be 229 //exclude edge pads
7390f655 230 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
231 Float_t valTmean=rocData->GetValue(irow,ipad);
892226be 232 //exclude values that are exactly 0
8166d5d7 233 if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
7390f655 234 rocOut->SetValue(irow,ipad,1);
892226be 235 ++noutliersCE;
236 }
237 // exclude channels with too large variations
238 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
7390f655 239 rocOut->SetValue(irow,ipad,1);
892226be 240 ++noutliersCE;
241 }
242 }
243 }
244 }
245 //perform fit
246 TMatrixD dummy;
2cb269df 247 Float_t chi2Af,chi2Cf;
248 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
249 chi2A=chi2Af;
250 chi2C=chi2Cf;
7390f655 251 if (!outCE) delete out;
892226be 252}
253//_____________________________________________________________________________________
254void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
255 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
256 Float_t &driftTimeA, Float_t &driftTimeC )
257{
258 //
259 // Calculate statistical information from the CE graphs for drift time and charge
260 //
261
262 //reset arrays
263 vecTEntries.ResizeTo(72);
264 vecTMean.ResizeTo(72);
265 vecTRMS.ResizeTo(72);
266 vecTMedian.ResizeTo(72);
267 vecQEntries.ResizeTo(72);
268 vecQMean.ResizeTo(72);
269 vecQRMS.ResizeTo(72);
270 vecQMedian.ResizeTo(72);
271 vecTEntries.Zero();
272 vecTMean.Zero();
273 vecTRMS.Zero();
274 vecTMedian.Zero();
275 vecQEntries.Zero();
276 vecQMean.Zero();
277 vecQRMS.Zero();
278 vecQMedian.Zero();
279 driftTimeA=0;
280 driftTimeC=0;
281 TObjArray *arrT=fCalibDB->GetCErocTtime();
282 TObjArray *arrQ=fCalibDB->GetCErocQtime();
283 if (arrT){
284 for (Int_t isec=0;isec<74;++isec){
285 TGraph *gr=(TGraph*)arrT->At(isec);
286 if (!gr) continue;
287 TVectorD values;
288 Int_t npoints = gr->GetN();
289 values.ResizeTo(npoints);
290 Int_t nused =0;
6e7d7dc4 291 //skip first points, theres always some problems with finding the CE position
292 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
293 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
892226be 294 values[nused]=gr->GetY()[ipoint];
295 nused++;
296 }
297 }
298 //
299 if (isec<72) vecTEntries[isec]= nused;
300 if (nused>1){
301 if (isec<72){
302 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
303 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
304 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
305 } else if (isec==72){
306 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
307 } else if (isec==73){
308 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
309 }
310 }
311 }
312 }
313 if (arrQ){
6e7d7dc4 314 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
892226be 315 TGraph *gr=(TGraph*)arrQ->At(isec);
316 if (!gr) continue;
317 TVectorD values;
318 Int_t npoints = gr->GetN();
319 values.ResizeTo(npoints);
320 Int_t nused =0;
321 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
817766d5 322 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
892226be 323 values[nused]=gr->GetY()[ipoint];
324 nused++;
325 }
326 }
327 //
6e7d7dc4 328 vecQEntries[isec]= nused;
892226be 329 if (nused>1){
330 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
331 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
332 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
333 }
334 }
335 }
336}
337
338//_____________________________________________________________________________________
339void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
340 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
93425263 341 Int_t &nonMaskedZero, Int_t &nNaN)
892226be 342{
343 //
344 // process noise data
345 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
346 // OROCs small pads [2] and OROCs large pads [3]
347 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
348 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
349 //
350
351 //set proper size and reset
352 const UInt_t infoSize=4;
353 vNoiseMean.ResizeTo(infoSize);
354 vNoiseMeanSenRegions.ResizeTo(infoSize);
355 vNoiseRMS.ResizeTo(infoSize);
356 vNoiseRMSSenRegions.ResizeTo(infoSize);
357 vNoiseMean.Zero();
358 vNoiseMeanSenRegions.Zero();
359 vNoiseRMS.Zero();
360 vNoiseRMSSenRegions.Zero();
361 nonMaskedZero=0;
93425263 362 nNaN=0;
892226be 363 //counters
364 TVectorD c(infoSize);
365 TVectorD cs(infoSize);
366 //tpc parameters
367 AliTPCParam par;
368 par.Update();
369 //retrieve noise and ALTRO data
370 if (!fPadNoise) return;
371 AliTPCCalROC *rocMasked=0x0;
372 //create IROC, OROC1, OROC2 and sensitive region masks
373 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
374 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
375 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
376 UInt_t nrows=noiseROC->GetNrows();
377 for (UInt_t irow=0;irow<nrows;++irow){
378 UInt_t npads=noiseROC->GetNPads(irow);
379 for (UInt_t ipad=0;ipad<npads;++ipad){
380 //don't use masked channels;
381 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
382 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
383 //check if noise==0
8166d5d7 384 if (noiseVal<kAlmost0) {
892226be 385 ++nonMaskedZero;
386 continue;
387 }
388 //check for nan
389 if ( !(noiseVal<10000000) ){
ef1807ef 390 AliInfo(Form("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal));
93425263 391 ++nNaN;
892226be 392 continue;
393 }
394 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
395 Int_t masksen=1; // sensitive pards are not masked (0)
396 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
397 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
398 //IROCs
399 if (irow>19&&irow<46){
400 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
401 }
402 Int_t type=1;
403 vNoiseMean[type]+=noiseVal;
404 vNoiseRMS[type]+=noiseVal*noiseVal;
405 ++c[type];
406 if (!masksen){
407 vNoiseMeanSenRegions[type]+=noiseVal;
408 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
409 ++cs[type];
410 }
411 } else {
412 //OROCs
413 //define sensive regions
414 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
415 if ( irow>75 ){
416 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
417 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
418 }
419 if ((Int_t)irow<par.GetNRowUp1()){
420 //OROC1
421 Int_t type=2;
422 vNoiseMean[type]+=noiseVal;
423 vNoiseRMS[type]+=noiseVal*noiseVal;
424 ++c[type];
425 if (!masksen){
426 vNoiseMeanSenRegions[type]+=noiseVal;
427 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
428 ++cs[type];
429 }
430 }else{
431 //OROC2
432 Int_t type=3;
433 vNoiseMean[type]+=noiseVal;
434 vNoiseRMS[type]+=noiseVal*noiseVal;
435 ++c[type];
436 if (!masksen){
437 vNoiseMeanSenRegions[type]+=noiseVal;
438 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
439 ++cs[type];
440 }
441 }
442 }
443 //whole tpc
444 Int_t type=0;
445 vNoiseMean[type]+=noiseVal;
446 vNoiseRMS[type]+=noiseVal*noiseVal;
447 ++c[type];
448 if (!masksen){
449 vNoiseMeanSenRegions[type]+=noiseVal;
450 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
451 ++cs[type];
452 }
453 }//end loop pads
454 }//end loop rows
455 }//end loop sectors (rocs)
456
457 //calculate mean and RMS
458 const Double_t verySmall=0.0000000001;
459 for (UInt_t i=0;i<infoSize;++i){
460 Double_t mean=0;
461 Double_t rms=0;
462 Double_t meanSen=0;
463 Double_t rmsSen=0;
464
465 if (c[i]>verySmall){
ef1807ef 466 AliInfo(Form("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]));
892226be 467 mean=vNoiseMean[i]/c[i];
468 rms=vNoiseRMS[i];
469 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
470 }
471 vNoiseMean[i]=mean;
472 vNoiseRMS[i]=rms;
473
474 if (cs[i]>verySmall){
475 meanSen=vNoiseMeanSenRegions[i]/cs[i];
476 rmsSen=vNoiseRMSSenRegions[i];
477 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
478 }
479 vNoiseMeanSenRegions[i]=meanSen;
480 vNoiseRMSSenRegions[i]=rmsSen;
481 }
482}
483
484//_____________________________________________________________________________________
abb20887 485void AliTPCcalibDButil::ProcessQAData(TVectorD &vQaOcc, TVectorD &vQaQtot,
486 TVectorD &vQaQmax)
487{
488 //
489 // process QA data
490 //
491 // vQaOcc/Qtot/Qmax contains the Mean occupancy/Qtot/Qmax for each sector
492 //
493
494
495 const UInt_t infoSize = 72;
496 //reset counters to error number
497 vQaOcc.ResizeTo(infoSize);
498 vQaOcc.Zero();
499 vQaQtot.ResizeTo(infoSize);
500 vQaQtot.Zero();
501 vQaQmax.ResizeTo(infoSize);
502 vQaQmax.Zero();
503 //counter
504 //retrieve pulser and ALTRO data
505
506 if (!fDataQA) {
507
508 AliInfo("No QA data");
509 return;
510 }
511 if (fDataQA->GetEventCounter()<=0) {
512
513 AliInfo("No QA data");
514 return; // no data processed
515 }
516 //
517 fDataQA->Analyse();
518
519 TVectorD normOcc(infoSize);
520 TVectorD normQ(infoSize);
521
522 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
523
ef1807ef 524 AliInfo(Form("Sector %d\n", isec));
abb20887 525 AliTPCCalROC* occupancyROC = fDataQA->GetNoThreshold()->GetCalROC(isec);
526 AliTPCCalROC* nclusterROC = fDataQA->GetNLocalMaxima()->GetCalROC(isec);
527 AliTPCCalROC* qROC = fDataQA->GetMeanCharge()->GetCalROC(isec);
528 AliTPCCalROC* qmaxROC = fDataQA->GetMaxCharge()->GetCalROC(isec);
529 if (!occupancyROC) continue;
530 if (!nclusterROC) continue;
531 if (!qROC) continue;
532 if (!qmaxROC) continue;
533
534 const UInt_t nchannels=occupancyROC->GetNchannels();
535
ef1807ef 536 AliInfo(Form("Nchannels %d\n", nchannels));
abb20887 537
538 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
539
540 vQaOcc[isec] += occupancyROC->GetValue(ichannel);
541 ++normOcc[isec];
542
543 Float_t nClusters = nclusterROC->GetValue(ichannel);
544 normQ[isec] += nClusters;
545 vQaQtot[isec]+=nClusters*qROC->GetValue(ichannel);
546 vQaQmax[isec]+=nClusters*qmaxROC->GetValue(ichannel);
547 }
548 }
549
550 //calculate mean values
551 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
552
553 if (normOcc[isec]>0) vQaOcc[isec] /= normOcc[isec];
554 else vQaOcc[isec] = 0;
555
556 if (normQ[isec]>0) {
557 vQaQtot[isec] /= normQ[isec];
558 vQaQmax[isec] /= normQ[isec];
559 }else {
560
561 vQaQtot[isec] = 0;
562 vQaQmax[isec] = 0;
563 }
564 }
565}
566
567//_____________________________________________________________________________________
892226be 568void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
569{
570 //
571 // Process the Pulser information
572 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
573 //
574
575 const UInt_t infoSize=4;
576 //reset counters to error number
577 vMeanTime.ResizeTo(infoSize);
578 vMeanTime.Zero();
579 //counter
580 TVectorD c(infoSize);
581 //retrieve pulser and ALTRO data
582 if (!fPulserTmean) return;
583 //
584 //get Outliers
585 AliTPCCalROC *rocOut=0x0;
586 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
587 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
588 if (!tmeanROC) continue;
589 rocOut=fPulserOutlier->GetCalROC(isec);
590 UInt_t nchannels=tmeanROC->GetNchannels();
591 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
592 if (rocOut && rocOut->GetValue(ichannel)) continue;
593 Float_t val=tmeanROC->GetValue(ichannel);
594 Int_t type=isec/18;
595 vMeanTime[type]+=val;
596 ++c[type];
597 }
598 }
599 //calculate mean
600 for (UInt_t itype=0; itype<infoSize; ++itype){
601 if (c[itype]>0) vMeanTime[itype]/=c[itype];
602 else vMeanTime[itype]=0;
603 }
604}
605//_____________________________________________________________________________________
606void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
607{
608 //
609 // Get Values from ALTRO configuration data
610 //
611 nMasked=-1;
612 if (!fALTROMasked) return;
613 nMasked=0;
614 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
615 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
616 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
617 if (rocMasked->GetValue(ichannel)) ++nMasked;
618 }
619 }
620}
621//_____________________________________________________________________________________
622void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
623{
624 //
625 // Proces Goofie values, return statistical information of the currently set goofieArray
626 // The meaning of the entries are given below
627 /*
628 1 TPC_ANODE_I_A00_STAT
629 2 TPC_DVM_CO2
630 3 TPC_DVM_DriftVelocity
631 4 TPC_DVM_FCageHV
632 5 TPC_DVM_GainFar
633 6 TPC_DVM_GainNear
634 7 TPC_DVM_N2
635 8 TPC_DVM_NumberOfSparks
636 9 TPC_DVM_PeakAreaFar
637 10 TPC_DVM_PeakAreaNear
638 11 TPC_DVM_PeakPosFar
639 12 TPC_DVM_PeakPosNear
640 13 TPC_DVM_PickupHV
641 14 TPC_DVM_Pressure
642 15 TPC_DVM_T1_Over_P
643 16 TPC_DVM_T2_Over_P
644 17 TPC_DVM_T_Over_P
645 18 TPC_DVM_TemperatureS1
646 */
647 if (!fGoofieArray){
648 Int_t nsensors=19;
649 vecEntries.ResizeTo(nsensors);
650 vecMedian.ResizeTo(nsensors);
651 vecMean.ResizeTo(nsensors);
652 vecRMS.ResizeTo(nsensors);
653 vecEntries.Zero();
654 vecMedian.Zero();
655 vecMean.Zero();
656 vecRMS.Zero();
657 return;
658 }
659 Double_t kEpsilon=0.0000000001;
660 Double_t kBig=100000000000.;
661 Int_t nsensors = fGoofieArray->NumSensors();
662 vecEntries.ResizeTo(nsensors);
663 vecMedian.ResizeTo(nsensors);
664 vecMean.ResizeTo(nsensors);
665 vecRMS.ResizeTo(nsensors);
666 TVectorF values;
667 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
668 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
669 if (gsensor && gsensor->GetGraph()){
670 Int_t npoints = gsensor->GetGraph()->GetN();
671 // filter zeroes
672 values.ResizeTo(npoints);
673 Int_t nused =0;
674 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
675 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
676 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
677 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
678 nused++;
679 }
680 }
681 //
682 vecEntries[isensor]= nused;
683 if (nused>1){
684 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
685 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
686 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
687 }
688 }
689 }
690}
7390f655 691//_____________________________________________________________________________________
692void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
693{
694 //
695 // check the variations of the pedestal data to the reference pedestal data
696 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
697 //
698 const Int_t npar=4;
699 TVectorF vThres(npar); //thresholds
700 Int_t nActive=0; //number of active channels
701
702 //reset and set thresholds
703 pedestalDeviations.ResizeTo(npar);
704 for (Int_t i=0;i<npar;++i){
705 pedestalDeviations.GetMatrixArray()[i]=0;
706 vThres.GetMatrixArray()[i]=(i+1)*.5;
707 }
708 //check all needed data is available
709 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
710 //loop over all channels
711 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
712 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
713 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
714 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
715 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
716 UInt_t nrows=mROC->GetNrows();
717 for (UInt_t irow=0;irow<nrows;++irow){
718 UInt_t npads=mROC->GetNPads(irow);
719 for (UInt_t ipad=0;ipad<npads;++ipad){
720 //don't use masked channels;
721 if (mROC ->GetValue(irow,ipad)) continue;
722 if (mRefROC->GetValue(irow,ipad)) continue;
723 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
724 for (Int_t i=0;i<npar;++i){
725 if (deviation>vThres[i])
726 ++pedestalDeviations.GetMatrixArray()[i];
727 }
728 ++nActive;
729 }//end ipad
730 }//ind irow
731 }//end isec
732 if (nActive>0){
733 for (Int_t i=0;i<npar;++i){
734 pedestalDeviations.GetMatrixArray()[i]/=nActive;
735 }
736 }
737}
738//_____________________________________________________________________________________
739void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
740{
741 //
742 // check the variations of the noise data to the reference noise data
743 // thresholds are 5, 10, 15 and 20 percent respectively.
744 //
745 const Int_t npar=4;
746 TVectorF vThres(npar); //thresholds
747 Int_t nActive=0; //number of active channels
748
749 //reset and set thresholds
750 noiseDeviations.ResizeTo(npar);
751 for (Int_t i=0;i<npar;++i){
752 noiseDeviations.GetMatrixArray()[i]=0;
753 vThres.GetMatrixArray()[i]=(i+1)*.05;
754 }
755 //check all needed data is available
756 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
757 //loop over all channels
758 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
759 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
760 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
761 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
762 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
763 UInt_t nrows=mROC->GetNrows();
764 for (UInt_t irow=0;irow<nrows;++irow){
765 UInt_t npads=mROC->GetNPads(irow);
766 for (UInt_t ipad=0;ipad<npads;++ipad){
767 //don't use masked channels;
768 if (mROC ->GetValue(irow,ipad)) continue;
769 if (mRefROC->GetValue(irow,ipad)) continue;
32100b1d 770 if (nRefROC->GetValue(irow,ipad)==0) continue;
7390f655 771 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
772 for (Int_t i=0;i<npar;++i){
773 if (deviation>vThres[i])
774 ++noiseDeviations.GetMatrixArray()[i];
775 }
776 ++nActive;
777 }//end ipad
778 }//ind irow
779 }//end isec
780 if (nActive>0){
781 for (Int_t i=0;i<npar;++i){
782 noiseDeviations.GetMatrixArray()[i]/=nActive;
783 }
784 }
785}
786//_____________________________________________________________________________________
787void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
788 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
789{
790 //
791 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
792 // thresholds are .5, 1, 5 and 10 percent respectively.
793 //
794 //
795 const Int_t npar=4;
796 TVectorF vThres(npar); //thresholds
797 Int_t nActive=0; //number of active channels
798
799 //reset and set thresholds
800 pulserQdeviations.ResizeTo(npar);
801 for (Int_t i=0;i<npar;++i){
802 pulserQdeviations.GetMatrixArray()[i]=0;
803 }
804 npadsOutOneTB=0;
805 npadsOffAdd=0;
806 varQMean=0;
807 vThres.GetMatrixArray()[0]=.005;
808 vThres.GetMatrixArray()[1]=.01;
809 vThres.GetMatrixArray()[2]=.05;
810 vThres.GetMatrixArray()[3]=.1;
811 //check all needed data is available
812 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
813 //
814 UpdateRefPulserOutlierMap();
815 //loop over all channels
732e90a8 816 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 817 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
818 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
819 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
732e90a8 820// AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
7390f655 821 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
822 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
823 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
d30aa177 824 Float_t ptmean=ptROC->GetMean(oROC);
7390f655 825 UInt_t nrows=mROC->GetNrows();
826 for (UInt_t irow=0;irow<nrows;++irow){
827 UInt_t npads=mROC->GetNPads(irow);
828 for (UInt_t ipad=0;ipad<npads;++ipad){
829 //don't use masked channels;
830 if (mROC ->GetValue(irow,ipad)) continue;
831 if (mRefROC->GetValue(irow,ipad)) continue;
832 //don't user edge pads
833 if (ipad==0||ipad==npads-1) continue;
834 //data
835 Float_t pq=pqROC->GetValue(irow,ipad);
836 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
837 Float_t pt=ptROC->GetValue(irow,ipad);
838// Float_t ptRef=ptRefROC->GetValue(irow,ipad);
839 //comparisons q
bb22b582 840 Float_t deviation=TMath::Abs(pqRef)>1e-20?TMath::Abs(pq/pqRef-1):-999;
7390f655 841 for (Int_t i=0;i<npar;++i){
842 if (deviation>vThres[i])
843 ++pulserQdeviations.GetMatrixArray()[i];
844 }
845 if (pqRef>11&&pq<11) ++npadsOffAdd;
846 varQMean+=pq-pqRef;
847 //comparisons t
d30aa177 848 if (TMath::Abs(pt-ptmean)>1) ++npadsOutOneTB;
7390f655 849 ++nActive;
850 }//end ipad
851 }//ind irow
852 }//end isec
853 if (nActive>0){
854 for (Int_t i=0;i<npar;++i){
855 pulserQdeviations.GetMatrixArray()[i]/=nActive;
856 varQMean/=nActive;
857 }
858 }
859}
892226be 860//_____________________________________________________________________________________
861void AliTPCcalibDButil::UpdatePulserOutlierMap()
862{
863 //
8166d5d7 864 // Update the outlier map of the pulser data
7390f655 865 //
866 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
867}
868//_____________________________________________________________________________________
869void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
870{
871 //
8166d5d7 872 // Update the outlier map of the pulser reference data
7390f655 873 //
874 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
875}
876//_____________________________________________________________________________________
877void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
878{
879 //
892226be 880 // Create a map that contains outliers from the Pulser calibration data.
881 // The outliers include masked channels, edge pads and pads with
882 // too large timing and charge variations.
7390f655 883 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
892226be 884 // those do not contain masked and edge pads
885 //
7390f655 886 if (!pulT||!pulQ) {
892226be 887 //reset map
7390f655 888 pulOut->Multiply(0.);
892226be 889 fNpulserOutliers=-1;
890 return;
891 }
892 AliTPCCalROC *rocMasked=0x0;
893 fNpulserOutliers=0;
894
895 //Create Outlier Map
896 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
7390f655 897 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
898 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
899 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
892226be 900 if (!tmeanROC||!qmeanROC) {
901 //reset outliers in this ROC
902 outROC->Multiply(0.);
903 continue;
904 }
905 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
906// Double_t dummy=0;
907// Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
908// Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
909 UInt_t nrows=tmeanROC->GetNrows();
910 for (UInt_t irow=0;irow<nrows;++irow){
911 UInt_t npads=tmeanROC->GetNPads(irow);
912 for (UInt_t ipad=0;ipad<npads;++ipad){
913 Int_t outlier=0,masked=0;
914 Float_t q=qmeanROC->GetValue(irow,ipad);
915 Float_t t=tmeanROC->GetValue(irow,ipad);
916 //masked channels are outliers
917 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
918 //edge pads are outliers
919 if (ipad==0||ipad==npads-1) masked=1;
920 //channels with too large charge or timing deviation from the meadian are outliers
921// if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
922 if (q<fPulQminLimit && !masked) outlier=1;
923 //check for nan
924 if ( !(q<10000000) || !(t<10000000)) outlier=1;
925 outROC->SetValue(irow,ipad,outlier+masked);
926 fNpulserOutliers+=outlier;
927 }
928 }
929 }
930}
931//_____________________________________________________________________________________
2cb269df 932AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
892226be 933{
934 //
935 // Create pad time0 object from pulser and/or CE data, depending on the selected model
936 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
937 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
732e90a8 938 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
892226be 939 //
2cb269df 940 // In case model 2 is invoked - gy arival time gradient is also returned
941 //
942 gyA=0;
943 gyC=0;
892226be 944 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
945 // decide between different models
946 if (model==0||model==1){
947 TVectorD vMean;
948 if (model==1) ProcessPulser(vMean);
949 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
950 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
951 if (!rocPulTmean) continue;
952 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
953 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
954 Float_t mean=rocPulTmean->GetMean(rocOut);
955 //treat case where a whole partition is masked
8166d5d7 956 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
892226be 957 if (model==1) {
958 Int_t type=isec/18;
959 mean=vMean[type];
960 }
961 UInt_t nrows=rocTime0->GetNrows();
962 for (UInt_t irow=0;irow<nrows;++irow){
963 UInt_t npads=rocTime0->GetNPads(irow);
964 for (UInt_t ipad=0;ipad<npads;++ipad){
965 Float_t time=rocPulTmean->GetValue(irow,ipad);
966 //in case of an outlier pad use the mean of the altro values.
967 //This should be the most precise guess in that case.
968 if (rocOut->GetValue(irow,ipad)) {
969 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
8166d5d7 970 if ( TMath::Abs(time)<kAlmost0 ) time=mean;
892226be 971 }
972 Float_t val=time-mean;
973 rocTime0->SetValue(irow,ipad,val);
974 }
975 }
976 }
2cb269df 977 } else if (model==2){
978 Double_t pgya,pgyc,pchi2a,pchi2c;
979 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
980 fCETmean->Add(padPulser,-1.);
732e90a8 981 TVectorD vA,vC;
982 AliTPCCalPad outCE("outCE","outCE");
983 Int_t nOut;
2cb269df 984 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
985 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
732e90a8 986// AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
2cb269df 987 if (!padFit) { delete padPulser; return 0;}
988 gyA=vA[2];
989 gyC=vC[2];
990 fCETmean->Add(padPulser,1.);
732e90a8 991 padTime0->Add(fCETmean);
2cb269df 992 padTime0->Add(padFit,-1);
993 delete padPulser;
732e90a8 994 TVectorD vFitROC;
995 TMatrixD mFitROC;
996 Float_t chi2;
997 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
998 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
999 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
1000 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
1001 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
1002 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
1003 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
1004 Float_t mean=rocPulTmean->GetMean(rocOutPul);
8166d5d7 1005 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
732e90a8 1006 UInt_t nrows=rocTime0->GetNrows();
1007 for (UInt_t irow=0;irow<nrows;++irow){
1008 UInt_t npads=rocTime0->GetNPads(irow);
1009 for (UInt_t ipad=0;ipad<npads;++ipad){
1010 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
1011 if (rocOutCE->GetValue(irow,ipad)){
1012 Float_t valOut=rocCEfit->GetValue(irow,ipad);
1013 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
1014 rocTime0->SetValue(irow,ipad,valOut);
1015 }
1016 }
1017 }
1018 delete rocCEfit;
1019 }
1020 delete padFit;
892226be 1021 }
2cb269df 1022 Double_t median = padTime0->GetMedian();
1023 padTime0->Add(-median); // normalize to median
892226be 1024 return padTime0;
1025}
1026//_____________________________________________________________________________________
98a4cc77 1027Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *const rocOut)
892226be 1028{
8166d5d7 1029 //
1030 // GetMeanAlto information
1031 //
892226be 1032 if (roc==0) return 0.;
1033 const Int_t sector=roc->GetSector();
1034 AliTPCROC *tpcRoc=AliTPCROC::Instance();
1035 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
1036 Float_t mean=0;
1037 Int_t n=0;
1038
1039 //loop over a small range around the requested pad (+-10 rows/pads)
1040 for (Int_t irow=row-10;irow<row+10;++irow){
1041 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
1042 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
1043 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
1044 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
1045 if (altroRoc!=altroCurr) continue;
1046 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
1047 Float_t val=roc->GetValue(irow,ipad);
1048 mean+=val;
1049 ++n;
1050 }
1051 }
1052 if (n>0) mean/=n;
1053 return mean;
1054}
7390f655 1055//_____________________________________________________________________________________
1056void AliTPCcalibDButil::SetRefFile(const char* filename)
1057{
1058 //
1059 // load cal pad objects form the reference file
1060 //
1061 TDirectory *currDir=gDirectory;
1062 TFile f(filename);
1063 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
1064 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
1065 //pulser data
1066 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
1067 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
1068 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
1069 //CE data
1070 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
1071 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
1072 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
1073 //Altro data
1074// fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
1075// fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
1076// fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
1077// fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
1078 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
1079 f.Close();
1080 currDir->cd();
1081}
949d8707 1082//_____________________________________________________________________________________
1083void AliTPCcalibDButil::UpdateRefDataFromOCDB()
1084{
1085 //
1086 // set reference data from OCDB Reference map
1087 //
1088 if (!fRefMap) {
1089 AliWarning("Referenc map not set!");
1090 return;
1091 }
1092
1093 TString cdbPath;
1094 AliCDBEntry* entry = 0x0;
1095 Bool_t hasAnyChanged=kFALSE;
1096
1097 //pedestals
1098 cdbPath="TPC/Calib/Pedestals";
1099 if (HasRefChanged(cdbPath.Data())){
1100 hasAnyChanged=kTRUE;
1101 //delete old entries
1102 if (fRefPedestals) delete fRefPedestals;
1103 if (fRefPedestalMasked) delete fRefPedestalMasked;
1104 fRefPedestals=fRefPedestalMasked=0x0;
1105 //get new entries
1106 entry=GetRefEntry(cdbPath.Data());
1107 if (entry){
1108 entry->SetOwner(kTRUE);
1109 fRefPedestals=GetRefCalPad(entry);
1110 delete entry;
1111 fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1112 }
1113 }
7390f655 1114
949d8707 1115 //noise
1116 cdbPath="TPC/Calib/PadNoise";
1117 if (HasRefChanged(cdbPath.Data())){
1118 hasAnyChanged=kTRUE;
1119 //delete old entry
1120 if (fRefPadNoise) delete fRefPadNoise;
1121 fRefPadNoise=0x0;
1122 //get new entry
1123 entry=GetRefEntry(cdbPath.Data());
1124 if (entry){
1125 entry->SetOwner(kTRUE);
1126 fRefPadNoise=GetRefCalPad(entry);
1127 delete entry;
1128 }
1129 }
1130
1131 //pulser
1132 cdbPath="TPC/Calib/Pulser";
1133 if (HasRefChanged(cdbPath.Data())){
1134 hasAnyChanged=kTRUE;
1135 //delete old entries
1136 if (fRefPulserTmean) delete fRefPulserTmean;
1137 if (fRefPulserTrms) delete fRefPulserTrms;
1138 if (fRefPulserQmean) delete fRefPulserQmean;
1139 if (fRefPulserMasked) delete fRefPulserMasked;
1140 fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1141 //get new entries
1142 entry=GetRefEntry(cdbPath.Data());
1143 if (entry){
1144 entry->SetOwner(kTRUE);
1145 fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1146 fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1147 fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1148 delete entry;
1149 fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1150 }
1151 }
2cb269df 1152
949d8707 1153 //ce
1154 cdbPath="TPC/Calib/CE";
1155 if (HasRefChanged(cdbPath.Data())){
1156 hasAnyChanged=kTRUE;
1157 //delete old entries
1158 if (fRefCETmean) delete fRefCETmean;
1159 if (fRefCETrms) delete fRefCETrms;
1160 if (fRefCEQmean) delete fRefCEQmean;
1161 if (fRefCEMasked) delete fRefCEMasked;
1162 fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1163 //get new entries
1164 entry=GetRefEntry(cdbPath.Data());
1165 if (entry){
1166 entry->SetOwner(kTRUE);
1167 fRefCETmean=GetRefCalPad(entry,"CETmean");
1168 fRefCETrms=GetRefCalPad(entry,"CETrms");
1169 fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1170 delete entry;
1171 fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1172 }
1173 }
1174
1175 //altro data
1176 cdbPath="TPC/Calib/AltroConfig";
1177 if (HasRefChanged(cdbPath.Data())){
1178 hasAnyChanged=kTRUE;
1179 //delete old entries
1180 if (fRefALTROFPED) delete fRefALTROFPED;
1181 if (fRefALTROZsThr) delete fRefALTROZsThr;
1182 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1183 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1184 if (fRefALTROMasked) delete fRefALTROMasked;
1185 fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1186 //get new entries
1187 entry=GetRefEntry(cdbPath.Data());
1188 if (entry){
1189 entry->SetOwner(kTRUE);
1190 fRefALTROFPED=GetRefCalPad(entry,"FPED");
1191 fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1192 fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1193 fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1194 fRefALTROMasked=GetRefCalPad(entry,"Masked");
1195 delete entry;
1196 }
1197 }
1198
1199 //raw data
1200 /*
1201 cdbPath="TPC/Calib/Raw";
1202 if (HasRefChanged(cdbPath.Data())){
1203 hasAnyChanged=kTRUE;
1204 //delete old entry
1205 if (fRefCalibRaw) delete fRefCalibRaw;
1206 //get new entry
1207 entry=GetRefEntry(cdbPath.Data());
1208 if (entry){
1209 entry->SetOwner(kTRUE);
1210 TObjArray *arr=(TObjArray*)entry->GetObject();
1211 if (!arr){
a980538f 1212 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1213 } else {
1214 fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1215 }
1216 }
1217 }
1218 */
2cb269df 1219
949d8707 1220 //data qa
1221 cdbPath="TPC/Calib/QA";
1222 if (HasRefChanged(cdbPath.Data())){
1223 hasAnyChanged=kTRUE;
1224 //delete old entry
1225 if (fRefDataQA) delete fRefDataQA;
1226 //get new entry
1227 entry=GetRefEntry(cdbPath.Data());
1228 if (entry){
1229 entry->SetOwner(kTRUE);
abb20887 1230 fRefDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1231 if (!fRefDataQA){
a980538f 1232 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1233 } else {
668e6ebb 1234 fRefDataQA=(AliTPCdataQA*)fRefDataQA->Clone();
949d8707 1235 }
1236 delete entry;
1237 }
1238 }
1239
1240
1241//update current reference maps
1242 if (hasAnyChanged){
1243 if (fCurrentRefMap) delete fCurrentRefMap;
1244 fCurrentRefMap=(TMap*)fRefMap->Clone();
1245 }
1246}
1247//_____________________________________________________________________________________
1248AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1249{
1250 //
1251 // TObjArray object type case
1252 // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1253 //
1254 AliTPCCalPad *pad=0x0;
1255 TObjArray *arr=(TObjArray*)entry->GetObject();
1256 if (!arr){
a980538f 1257 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1258 return pad;
1259 }
1260 pad=(AliTPCCalPad*)arr->FindObject(objName);
1261 if (!pad) {
a980538f 1262 AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
949d8707 1263 return pad;
1264 }
1265 return (AliTPCCalPad*)pad->Clone();
1266}
1267//_____________________________________________________________________________________
1268AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1269{
1270 //
1271 // AliTPCCalPad object type case
1272 // cast object to a calPad and store it in 'pad'
1273 //
1274 AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1275 if (!pad) {
a980538f 1276 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
949d8707 1277 return 0x0;
1278 }
1279 pad=(AliTPCCalPad*)pad->Clone();
1280 return pad;
1281}
1282//_____________________________________________________________________________________
1283AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1284{
1285 //
1286 // set altro masked channel map for 'cdbPath'
1287 //
1288 AliTPCCalPad* pad=0x0;
1289 const Int_t run=GetReferenceRun(cdbPath);
1290 if (run<0) {
a980538f 1291 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1292 return pad;
1293 }
1294 AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1295 if (!entry) {
a980538f 1296 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1297 return pad;
1298 }
1299 pad=GetRefCalPad(entry,"Masked");
1300 if (pad) pad->SetNameTitle(name,name);
1301 entry->SetOwner(kTRUE);
1302 delete entry;
1303 return pad;
1304}
1305//_____________________________________________________________________________________
1306void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1307 //
1308 // Get Reference map
1309 //
1310 if (run<0) run=fCalibDB->GetRun();
1311 TString cdbPath="TPC/Calib/Ref";
1312 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1313 if (!entry) {
a980538f 1314 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
949d8707 1315 fRefMap=0;
1316 return;
1317 }
1318 entry->SetOwner(kTRUE);
1319 fRefMap=(TMap*)(entry->GetObject());
1320 AliCDBId &id=entry->GetId();
1321 fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1322}
1323//_____________________________________________________________________________________
1324Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1325{
1326 //
1327 // check whether a reference cdb entry has changed
1328 //
1329 if (!fCurrentRefMap) return kTRUE;
1330 if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1331 return kFALSE;
1332}
1333//_____________________________________________________________________________________
1334AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1335{
1336 //
1337 // get the reference AliCDBEntry for 'cdbPath'
1338 //
1339 const Int_t run=GetReferenceRun(cdbPath);
1340 if (run<0) {
a980538f 1341 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1342 return 0;
1343 }
1344 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1345 if (!entry) {
a980538f 1346 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
949d8707 1347 return 0;
1348 }
1349 return entry;
1350}
1351//_____________________________________________________________________________________
26d8859c 1352Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
949d8707 1353 //
1354 // Get reference run number for the specified OCDB path
1355 //
1356 if (!fCurrentRefMap) return -2;
1357 TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1358 if (!str) return -2;
58affc42 1359 return (Int_t)str->GetString().Atoi();
949d8707 1360}
1361//_____________________________________________________________________________________
26d8859c 1362Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
949d8707 1363 //
1364 // Get reference run number for the specified OCDB path
1365 //
1366 if (!fRefMap) return -1;
1367 TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1368 if (!str) return -1;
58affc42 1369 return (Int_t)str->GetString().Atoi();
949d8707 1370}
1371//_____________________________________________________________________________________
8166d5d7 1372AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
2cb269df 1373 //
1374 // Author: marian.ivanov@cern.ch
1375 //
1376 // Create outlier map for CE study
1377 // Parameters:
1378 // Return value - outlyer map
1379 // noutlyersCE - number of outlyers
1380 // minSignal - minimal total Q signal
1381 // cutRMSMin - minimal width of the signal in respect to the median
1382 // cutRMSMax - maximal width of the signal in respect to the median
1383 // cutMaxDistT - maximal deviation from time median per chamber
1384 //
1385 // Outlyers criteria:
1386 // 0. Exclude masked pads
1387 // 1. Exclude first two rows in IROC and last two rows in OROC
1388 // 2. Exclude edge pads
1389 // 3. Exclude channels with too large variations
1390 // 4. Exclude pads with too small signal
1391 // 5. Exclude signal with outlyers RMS
1392 // 6. Exclude channels to far from the chamber median
1393 noutliersCE=0;
1394 //create outlier map
1395 AliTPCCalPad *out=ceOut;
1396 if (!out) out= new AliTPCCalPad("outCE","outCE");
1397 AliTPCCalROC *rocMasked=0x0;
1398 if (!fCETmean) return 0;
1399 if (!fCETrms) return 0;
1400 if (!fCEQmean) return 0;
1401 //
1402 //loop over all channels
1403 //
1404 Double_t rmsMedian = fCETrms->GetMedian();
1405 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1406 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
9f98a33d 1407 if (!rocData) continue;
2cb269df 1408 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1409 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1410 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
1411 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1412 Double_t trocMedian = rocData->GetMedian();
1413 //
768e46f1 1414 if (!rocData || !rocCEQ || !rocCETrms || !rocData) {
2cb269df 1415 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1416 rocOut->Add(1.);
1417 continue;
1418 }
1419 //
1420 //select outliers
1421 UInt_t nrows=rocData->GetNrows();
1422 for (UInt_t irow=0;irow<nrows;++irow){
1423 UInt_t npads=rocData->GetNPads(irow);
1424 for (UInt_t ipad=0;ipad<npads;++ipad){
1425 rocOut->SetValue(irow,ipad,0);
1426 Float_t valTmean=rocData->GetValue(irow,ipad);
1427 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1428 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1429 //0. exclude masked pads
1430 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1431 rocOut->SetValue(irow,ipad,1);
1432 continue;
1433 }
1434 //1. exclude first two rows in IROC and last two rows in OROC
1435 if (iroc<36){
1436 if (irow<2) rocOut->SetValue(irow,ipad,1);
1437 } else {
1438 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1439 }
1440 //2. exclude edge pads
1441 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1442 //exclude values that are exactly 0
8166d5d7 1443 if ( TMath::Abs(valTmean)<kAlmost0) {
2cb269df 1444 rocOut->SetValue(irow,ipad,1);
1445 ++noutliersCE;
1446 }
1447 //3. exclude channels with too large variations
1448 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1449 rocOut->SetValue(irow,ipad,1);
1450 ++noutliersCE;
1451 }
1452 //
1453 //4. exclude channels with too small signal
1454 if (valQmean<minSignal) {
1455 rocOut->SetValue(irow,ipad,1);
1456 ++noutliersCE;
1457 }
1458 //
1459 //5. exclude channels with too small rms
1460 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1461 rocOut->SetValue(irow,ipad,1);
1462 ++noutliersCE;
1463 }
1464 //
1465 //6. exclude channels to far from the chamber median
1466 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1467 rocOut->SetValue(irow,ipad,1);
1468 ++noutliersCE;
1469 }
1470 }
1471 }
1472 }
1473 //
1474 return out;
1475}
1476
1477
8166d5d7 1478AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
2cb269df 1479 //
1480 // Author: marian.ivanov@cern.ch
1481 //
1482 // Create outlier map for Pulser
1483 // Parameters:
1484 // Return value - outlyer map
1485 // noutlyersPulser - number of outlyers
1486 // cutTime - absolute cut - distance to the median of chamber
1487 // cutnRMSQ - nsigma cut from median q distribution per chamber
1488 // cutnRMSrms - nsigma cut from median rms distribution
1489 // Outlyers criteria:
1490 // 0. Exclude masked pads
1491 // 1. Exclude time outlyers (default 3 time bins)
1492 // 2. Exclude q outlyers (default 5 sigma)
1493 // 3. Exclude rms outlyers (default 5 sigma)
1494 noutliersPulser=0;
1495 AliTPCCalPad *out=pulserOut;
1496 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1497 AliTPCCalROC *rocMasked=0x0;
1498 if (!fPulserTmean) return 0;
1499 if (!fPulserTrms) return 0;
1500 if (!fPulserQmean) return 0;
1501 //
1502 //loop over all channels
1503 //
1504 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1505 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1506 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1507 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1508 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1509 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1510 //
1511 Double_t rocMedianT = rocData->GetMedian();
1512 Double_t rocMedianQ = rocPulserQ->GetMedian();
1513 Double_t rocRMSQ = rocPulserQ->GetRMS();
1514 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1515 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1516 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1517 rocOut->SetValue(ichannel,0);
1518 Float_t valTmean=rocData->GetValue(ichannel);
1519 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1520 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
f03e487b 1521 Float_t valMasked =0;
1522 if (rocMasked) valMasked = rocMasked->GetValue(ichannel);
2cb269df 1523 Int_t isOut=0;
f03e487b 1524 if (valMasked>0.5) isOut=1;
2cb269df 1525 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1526 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1527 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1528 rocOut->SetValue(ichannel,isOut);
1529 if (isOut) noutliersPulser++;
1530 }
1531 }
1532 return out;
1533}
1534
1535
1536AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1537 //
1538 // Author : Marian Ivanov
1539 // Create pad time0 correction map using information from the CE and from pulser
1540 //
1541 //
1542 // Return PadTime0 to be used for time0 relative alignment
1543 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1544 // using $ALICE_ROOT/TPC/script/gui application
1545 //
1546 // fitResultsA - fitParameters A side
1547 // fitResultsC - fitParameters C side
1548 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1549 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1550 //
1551 //
1552 // Algorithm:
1553 // 1. Find outlier map for CE
1554 // 2. Find outlier map for Pulser
1555 // 3. Replace outlier by median at given sector (median without outliers)
1556 // 4. Substract from the CE data pulser
1557 // 5. Fit the CE with formula
1558 // 5.1) (IROC-OROC) offset
1559 // 5.2) gx
1560 // 5.3) gy
1561 // 5.4) (lx-xmid)
1562 // 5.5) (IROC-OROC)*(lx-xmid)
1563 // 5.6) (ly/lx)^2
1564 // 6. Substract gy fit dependence from the CE data
1565 // 7. Add pulser back to CE data
1566 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1567 // 9. return CE data
1568 //
1569 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1570 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1571
1572 // fit formula
1573 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)";
1574 // output for fit formula
1575 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)";
1576 // gy part of formula
1577 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)";
1578 //
1579 //
1580 if (!fCETmean) return 0;
1581 Double_t pgya,pgyc,pchi2a,pchi2c;
1582 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1583 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1584
1585 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1586 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1587 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1588 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1589 padPulser->SetName("padPulser");
1590 padPulserOut->SetName("padPulserOut");
1591 padCE->SetName("padCE");
1592 padCEIn->SetName("padCEIn");
1593 padCEOut->SetName("padCEOut");
1594 padOut->SetName("padOut");
1595
1596 //
1597 // make combined outlyers map
1598 // and replace outlyers in maps with median for chamber
1599 //
1600 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1601 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1602 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1603 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1604 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1605 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1606 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1607 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1608 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1609 if (rocPulserOut->GetValue(ichannel)>0) {
1610 rocPulser->SetValue(ichannel,pulserMedian);
1611 rocOut->SetValue(ichannel,1);
1612 }
1613 if (rocCEOut->GetValue(ichannel)>0) {
1614 rocCE->SetValue(ichannel,ceMedian);
1615 rocOut->SetValue(ichannel,1);
1616 }
1617 }
1618 }
1619 //
1620 // remove pulser time 0
1621 //
1622 padCE->Add(padPulser,-1);
1623 //
1624 // Make fits
1625 //
1626 TMatrixD dummy;
1627 Float_t chi2Af,chi2Cf;
1628 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1629 chi2A=chi2Af;
1630 chi2C=chi2Cf;
1631 //
1632 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1633 padCEFitGY->SetName("padCEFitGy");
1634 //
1635 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1636 padCEFit->SetName("padCEFit");
1637 //
1638 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1639 padCEDiff->SetName("padCEDiff");
1640 padCEDiff->Add(padCEFit,-1.);
1641 //
1642 //
1643 padCE->Add(padCEFitGY,-1.);
1644
1645 padCE->Add(padPulser,1.);
1646 Double_t padmedian = padCE->GetMedian();
1647 padCE->Add(-padmedian); // normalize to median
1648 //
1649 // Replace outliers by fit value - median of diff per given chamber -GY fit
1650 //
1651 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1652 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1653 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1654 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1655 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1656 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1657 //
1658 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1659 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1660 if (rocOut->GetValue(ichannel)==0) continue;
1661 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1662 rocCE->SetValue(ichannel,value);
1663 }
1664 }
1665 //
1666 //
1667 if (dumpfile){
1668 //dump to the file - result can be visualized
1669 AliTPCPreprocessorOnline preprocesor;
1670 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1671 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1672 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1673 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1674 //
1675 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1676 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1677 //
1678 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1679 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1680 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1681 preprocesor.DumpToFile(dumpfile);
1682 }
1683 delete padPulser;
1684 delete padPulserOut;
1685 delete padCEIn;
1686 delete padCEOut;
1687 delete padOut;
1688 delete padCEDiff;
1689 delete padCEFitGY;
1690 return padCE;
1691}
1692
817766d5 1693
1694
1695
1696
1697Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1698 //
1699 // find the closest point to xref in x direction
1700 // return dx and value
5647625c 1701 dx = 0;
1702 y = 0;
1703
1704 if(!graph) return 0;
1705 if(graph->GetN() < 1) return 0;
1706
817766d5 1707 Int_t index=0;
1708 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1709 if (index<0) index=0;
5647625c 1710 if(graph->GetN()==1) {
1711 dx = xref-graph->GetX()[index];
1712 }
1713 else {
1714 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1715 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1716 dx = xref-graph->GetX()[index];
1717 }
817766d5 1718 y = graph->GetY()[index];
1719 return index;
1720}
1721
817766d5 1722Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1723 //
1724 // Get the correction of the trigger offset
1725 // combining information from the laser track calibration
1726 // and from cosmic calibration
1727 //
1728 // run - run number
1729 // timeStamp - tim stamp in seconds
1730 // deltaT - integration period to calculate offset
1731 // deltaTLaser -max validity of laser data
1732 // valType - 0 - median, 1- mean
1733 //
1734 // Integration vaues are just recomendation - if not possible to get points
1735 // automatically increase the validity by factor 2
1736 // (recursive algorithm until one month of data taking)
1737 //
1738 //
1739 const Float_t kLaserCut=0.0005;
aec07590 1740 const Int_t kMaxPeriod=3600*24*30*12; // one year max
817766d5 1741 const Int_t kMinPoints=20;
1742 //
1743 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1744 if (!array) {
1745 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1746 }
1747 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1748 if (!array) return 0;
1749 //
1750 TGraphErrors *laserA[3]={0,0,0};
1751 TGraphErrors *laserC[3]={0,0,0};
1752 TGraphErrors *cosmicAll=0;
1753 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1754 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1755 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1756 //
1757 //
1758 if (!cosmicAll) return 0;
1759 Int_t nmeasC=cosmicAll->GetN();
1760 Float_t *tdelta = new Float_t[nmeasC];
1761 Int_t nused=0;
1762 for (Int_t i=0;i<nmeasC;i++){
1763 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1764 Float_t ccosmic=cosmicAll->GetY()[i];
1765 Double_t yA=0,yC=0,dA=0,dC=0;
1766 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1767 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1768 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1769 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1770 //
1771 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1772 Float_t claser=0;
1773 if (TMath::Abs(yA-yC)<kLaserCut) {
1774 claser=(yA-yC)*0.5;
1775 }else{
1776 if (i%2==0) claser=yA;
1777 if (i%2==1) claser=yC;
1778 }
1779 tdelta[nused]=ccosmic-claser;
1780 nused++;
1781 }
768e46f1 1782 if (nused<kMinPoints &&deltaT<kMaxPeriod) {
1783 delete [] tdelta;
1784 return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1785 }
aec07590 1786 if (nused<kMinPoints) {
f03e487b 1787 delete [] tdelta;
6157d837 1788 //AliWarning("AliFatal: No time offset calibration available\n");
aec07590 1789 return 0;
1790 }
817766d5 1791 Double_t median = TMath::Median(nused,tdelta);
1792 Double_t mean = TMath::Mean(nused,tdelta);
a9f487a1 1793 delete [] tdelta;
817766d5 1794 return (valType==0) ? median:mean;
1795}
1796
a23ba1c3 1797Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
817766d5 1798 //
1799 // Get the correction of the drift velocity
1800 // combining information from the laser track calibration
1801 // and from cosmic calibration
1802 //
a23ba1c3 1803 // dist - return value - distance to closest point in graph
817766d5 1804 // run - run number
1805 // timeStamp - tim stamp in seconds
1806 // deltaT - integration period to calculate time0 offset
1807 // deltaTLaser -max validity of laser data
1808 // valType - 0 - median, 1- mean
1809 //
1810 // Integration vaues are just recomendation - if not possible to get points
1811 // automatically increase the validity by factor 2
1812 // (recursive algorithm until one month of data taking)
1813 //
1814 //
1815 //
1816 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1817 if (!array) {
1818 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1819 }
1820 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1821 if (!array) return 0;
1822 TGraphErrors *cosmicAll=0;
1823 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1824 if (!cosmicAll) return 0;
a23ba1c3 1825 Double_t grY=0;
1826 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1827
817766d5 1828 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
5647625c 1829 Double_t vcosmic = AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1e722a63 1830 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
817766d5 1831 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
cc65e4f5 1832 return vcosmic-t0;
817766d5 1833
1834 /*
1835 Example usage:
1836
1837 Int_t run=89000
1838 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1839 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1840 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1841 //
1842 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1843 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1844 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1845 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1846
1847 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1848 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1849
1850 */
1851
1852}
1853
949d8707 1854const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1855{
1856 //
1857 // Create a default name for the gui file
1858 //
1859
1860 return Form("guiRefTreeRun%s.root",GetRefValidity());
1861}
1862
1863Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1864{
1865 //
1866 // Create a gui reference tree
1867 // if dirname and filename are empty default values will be used
1868 // this is the recommended way of using this function
1869 // it allows to check whether a file with the given run validity alredy exists
1870 //
1871 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1872 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1873 return kFALSE;
1874 }
1875
1876 TString file=filename;
1877 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1878
1879 AliTPCPreprocessorOnline prep;
1880 //noise and pedestals
1881 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1882 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1883 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1884 //pulser data
1885 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1886 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1887 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1888 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1889 //CE data
1890 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1891 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1892 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1893 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1894 //Altro data
1895 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1896 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1897 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1898 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1899 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1900 //QA
1901 AliTPCdataQA *dataQA=fRefDataQA;
1902 if (dataQA) {
1903 if (dataQA->GetNLocalMaxima())
1904 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1905 if (dataQA->GetMaxCharge())
1906 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1907 if (dataQA->GetMeanCharge())
1908 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1909 if (dataQA->GetNoThreshold())
1910 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1911 if (dataQA->GetNTimeBins())
1912 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1913 if (dataQA->GetNPads())
1914 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1915 if (dataQA->GetTimePosition())
1916 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1917 }
1918 prep.DumpToFile(file.Data());
1919 return kTRUE;
1920}
1921
a23ba1c3 1922Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1923 //
78f17711 1924 // Get the correction of the drift velocity using the offline laser tracks calbration
a23ba1c3 1925 //
1926 // run - run number
1927 // timeStamp - tim stamp in seconds
1928 // deltaT - integration period to calculate time0 offset
1929 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1930 // Note in case no data form both A and C side - the value from active side used
1931 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
78f17711 1932
1933 return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1934}
1935
1936Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1937 //
1938 // Get the correction of the drift velocity using the online laser tracks calbration
1939 //
1940 // run - run number
1941 // timeStamp - tim stamp in seconds
1942 // deltaT - integration period to calculate time0 offset
1943 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1944 // Note in case no data form both A and C side - the value from active side used
1945 TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1946
1947 Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1948 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1949 if (!param) return 0;
1950
1951 //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1952 dv*=param->GetDriftV()/2.61301900000000000e+06;
1953 if (dv>1e-20) dv=1/dv-1;
1954 else return 0;
1955 // T/P correction
1956 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1957
1958 AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
2cd3a697 1959 AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
78f17711 1960
1961 Double_t corrPTA=0;
1962 Double_t corrPTC=0;
1963
1964 if (temp&&press) {
1965 AliTPCCalibVdrift corr(temp,press,0);
1966 corrPTA=corr.GetPTRelative(timeStamp,0);
1967 corrPTC=corr.GetPTRelative(timeStamp,1);
1968 }
1969
1970 if (side==0) dv -= corrPTA;
1971 if (side==1) dv -= corrPTC;
1972 if (side==2) dv -= (corrPTA+corrPTC)/2;
1973
1974 return dv;
1975}
1976
1977Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1978 Int_t side, TObjArray * const array){
1979 //
1980 // common drift velocity retrieval for online and offline method
1981 //
a23ba1c3 1982 TGraphErrors *grlaserA=0;
1983 TGraphErrors *grlaserC=0;
1984 Double_t vlaserA=0, vlaserC=0;
1985 if (!array) return 0;
1986 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1987 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1988 Double_t deltaY;
5647625c 1989 if (grlaserA && grlaserA->GetN()>0) {
a23ba1c3 1990 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1991 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1992 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1993 }
5647625c 1994 if (grlaserC && grlaserC->GetN()>0) {
a23ba1c3 1995 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1996 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1997 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1998 }
1999 if (side==0) return vlaserA;
2000 if (side==1) return vlaserC;
2001 Double_t mdrift=(vlaserA+vlaserC)*0.5;
2002 if (!grlaserA) return vlaserC;
2003 if (!grlaserC) return vlaserA;
2004 return mdrift;
2005}
2006
2007
2008Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
2009 //
2010 // Get the correction of the drift velocity using the CE laser data
2011 // combining information from the CE, laser track calibration
2012 // and P/T calibration
2013 //
2014 // run - run number
2015 // timeStamp - tim stamp in seconds
2016 // deltaT - integration period to calculate time0 offset
2017 // side - 0 - A side, 1 - C side, 2 - mean from both sides
2018 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
0d1b4cf8 2019 if (!arrT) return 0;
a23ba1c3 2020 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
2021 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2022 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2023 //
2024 //
2025 Double_t corrPTA = 0, corrPTC=0;
2026 Double_t ltime0A = 0, ltime0C=0;
2027 Double_t gry=0;
2028 Double_t corrA=0, corrC=0;
2029 Double_t timeA=0, timeC=0;
5647625c 2030 const Double_t kEpsilon = 0.00001;
a23ba1c3 2031 TGraph *graphA = (TGraph*)arrT->At(72);
2032 TGraph *graphC = (TGraph*)arrT->At(73);
2033 if (!graphA && !graphC) return 0.;
2034 if (graphA &&graphA->GetN()>0) {
2035 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2036 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2037 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2038 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2039 if(ltime0A < kEpsilon) return 0;
a23ba1c3 2040 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2041 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2042 corrA-=corrPTA;
2043 }
2044 if (graphC&&graphC->GetN()>0){
2045 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2046 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2047 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2048 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
5647625c 2049 if(ltime0C < kEpsilon) return 0;
2050if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1e722a63 2051 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
a23ba1c3 2052 corrC-=corrPTC;
2053 }
2054
2055 if (side ==0 ) return corrA;
2056 if (side ==1 ) return corrC;
2057 Double_t corrM= (corrA+corrC)*0.5;
2058 if (!graphA) corrM=corrC;
2059 if (!graphC) corrM=corrA;
2060 return corrM;
2061}
2062
cc65e4f5 2063Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2064 //
2065 // return drift velocity using the TPC-ITS matchin method
2066 // return also distance to the closest point
2067 //
2068 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2069 TGraphErrors *graph=0;
2070 dist=0;
2071 if (!array) return 0;
5647625c 2072 //array->ls();
cc65e4f5 2073 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
3760239f 2074 if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_TOFB_TPC_DRIFTVD");
cc65e4f5 2075 if (!graph) return 0;
2076 Double_t deltaY;
2077 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2078 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2079 return value;
2080}
2081
2082Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2083 //
2084 // Get time dependent time 0 (trigger delay in cm) correction
2085 // Arguments:
2086 // timestamp - timestamp
2087 // run - run number
2088 //
2089 // Notice - Extrapolation outside of calibration range - using constant function
2090 //
2091 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2092 TGraphErrors *graph=0;
2093 dist=0;
2094 if (!array) return 0;
3760239f 2095 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
2096 if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
cc65e4f5 2097 if (!graph) return 0;
2098 Double_t deltaY;
2099 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2100 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2101 return value;
2102}
2103
2104
a23ba1c3 2105
2106
2107
2108Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2109 //
2110 // VERY obscure method - we need something in framework
2111 // Find the TPC runs with temperature OCDB entry
2112 // cache the start and end of the run
2113 //
2114 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2115 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2116 if (!storage) return 0;
2117 TString path=storage->GetURI();
2118 TString runsT;
2119 {
2120 TString command;
2121 if (path.Contains("local")){ // find the list if local system
2122 path.ReplaceAll("local://","");
2123 path+="TPC/Calib/Temperature";
2124 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2125 }
2126 runsT=gSystem->GetFromPipe(command);
2127 }
2128 TObjArray *arr= runsT.Tokenize("r");
2129 if (!arr) return 0;
2130 //
2131 TArrayI indexes(arr->GetEntries());
2132 TArrayI runs(arr->GetEntries());
2133 Int_t naccept=0;
2134 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2135 Int_t irunN = atoi(arr->At(irun)->GetName());
2136 if (irunN<startRun) continue;
2137 if (irunN>stopRun) continue;
2138 runs[naccept]=irunN;
2139 naccept++;
2140 }}
09d5920f 2141 delete arr;
a23ba1c3 2142 fRuns.Set(naccept);
2143 fRunsStart.Set(fRuns.fN);
2144 fRunsStop.Set(fRuns.fN);
2145 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2146 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2147
2148 //
2149 AliCDBEntry * entry = 0;
2150 {for (Int_t irun=0;irun<fRuns.fN; irun++){
2151 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2152 if (!entry) continue;
2153 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2154 if (!tmpRun) continue;
2155 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2156 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
ef1807ef 2157 //AliInfo(Form("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec()));
a23ba1c3 2158 }}
2159 return fRuns.fN;
2160}
2161
2162
2163Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2164 //
2165 // binary search - find the run for given time stamp
2166 //
2167 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2168 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2169 Int_t cindex = -1;
2170 for (Int_t index=index0; index<=index1; index++){
2171 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2172 if (debug) {
ef1807ef 2173 AliInfo(Form("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime));
a23ba1c3 2174 }
2175 }
2176 if (cindex<0) cindex =(index0+index1)/2;
2177 if (cindex<0) {
2178 return 0;
2179 }
2180 return fRuns[cindex];
2181}
2182
2183
2184
2185
2186
2187TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2188 //
2189 // filter outlyer measurement
2190 // Only points around median +- sigmaCut filtered
2191 if (!graph) return 0;
2192 Int_t kMinPoints=2;
2193 Int_t npoints0 = graph->GetN();
2194 Int_t npoints=0;
2195 Float_t rmsY=0;
a23ba1c3 2196 //
2197 //
2198 if (npoints0<kMinPoints) return 0;
a9f487a1 2199
2200 Double_t *outx=new Double_t[npoints0];
2201 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2202 for (Int_t iter=0; iter<3; iter++){
2203 npoints=0;
2204 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2205 if (graph->GetY()[ipoint]==0) continue;
2206 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2207 outx[npoints] = graph->GetX()[ipoint];
2208 outy[npoints] = graph->GetY()[ipoint];
2209 npoints++;
2210 }
2211 if (npoints<=1) break;
2212 medianY =TMath::Median(npoints,outy);
2213 rmsY =TMath::RMS(npoints,outy);
2214 }
2215 TGraph *graphOut=0;
2216 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2217 delete [] outx;
2218 delete [] outy;
a23ba1c3 2219 return graphOut;
2220}
2221
2222
2223TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2224 //
2225 // filter outlyer measurement
2226 // Only points around median +- cut filtered
2227 if (!graph) return 0;
2228 Int_t kMinPoints=2;
2229 Int_t npoints0 = graph->GetN();
2230 Int_t npoints=0;
829455ad 2231// Float_t rmsY=0;
a23ba1c3 2232 //
2233 //
2234 if (npoints0<kMinPoints) return 0;
a9f487a1 2235
2236 Double_t *outx=new Double_t[npoints0];
2237 Double_t *outy=new Double_t[npoints0];
a23ba1c3 2238 for (Int_t iter=0; iter<3; iter++){
2239 npoints=0;
2240 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2241 if (graph->GetY()[ipoint]==0) continue;
2242 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2243 outx[npoints] = graph->GetX()[ipoint];
2244 outy[npoints] = graph->GetY()[ipoint];
2245 npoints++;
2246 }
2247 if (npoints<=1) break;
2248 medianY =TMath::Median(npoints,outy);
829455ad 2249// rmsY =TMath::RMS(npoints,outy);
a23ba1c3 2250 }
2251 TGraph *graphOut=0;
2252 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
a9f487a1 2253 delete [] outx;
2254 delete [] outy;
a23ba1c3 2255 return graphOut;
2256}
2257
2258
2259
8166d5d7 2260TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
a23ba1c3 2261 //
2262 // filter outlyer measurement
2263 // Only points with normalized errors median +- sigmaCut filtered
2264 //
2265 Int_t kMinPoints=10;
2266 Int_t npoints0 = graph->GetN();
2267 Int_t npoints=0;
2268 Float_t medianErr=0, rmsErr=0;
a9f487a1 2269 //
2270 //
2271 if (npoints0<kMinPoints) return 0;
2272
a23ba1c3 2273 Double_t *outx=new Double_t[npoints0];
2274 Double_t *outy=new Double_t[npoints0];
2275 Double_t *erry=new Double_t[npoints0];
2276 Double_t *nerry=new Double_t[npoints0];
2277 Double_t *errx=new Double_t[npoints0];
a9f487a1 2278
a23ba1c3 2279 for (Int_t iter=0; iter<3; iter++){
2280 npoints=0;
2281 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2282 nerry[npoints] = graph->GetErrorY(ipoint);
2283 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2284 erry[npoints] = graph->GetErrorY(ipoint);
2285 outx[npoints] = graph->GetX()[ipoint];
2286 outy[npoints] = graph->GetY()[ipoint];
2287 errx[npoints] = graph->GetErrorY(ipoint);
2288 npoints++;
2289 }
2290 if (npoints==0) break;
2291 medianErr=TMath::Median(npoints,erry);
2292 medianY =TMath::Median(npoints,outy);
2293 rmsErr =TMath::RMS(npoints,erry);
2294 }
2295 TGraphErrors *graphOut=0;
2296 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2297 delete []outx;
2298 delete []outy;
a23ba1c3 2299 delete []erry;
a9f487a1 2300 delete []nerry;
2301 delete []errx;
a23ba1c3 2302 return graphOut;
2303}
2304
2305
2306void AliTPCcalibDButil::Sort(TGraph *graph){
2307 //
2308 // sort array - neccessay for approx
2309 //
2310 Int_t npoints = graph->GetN();
2311 Int_t *indexes=new Int_t[npoints];
2312 Double_t *outx=new Double_t[npoints];
2313 Double_t *outy=new Double_t[npoints];
2314 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2315 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2316 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2317 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2318 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2319
a9f487a1 2320 delete [] indexes;
2321 delete [] outx;
2322 delete [] outy;
a23ba1c3 2323}
2324void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2325 //
2326 // smmoth graph - mean on the interval
2327 //
2328 Sort(graph);
2329 Int_t npoints = graph->GetN();
2330 Double_t *outy=new Double_t[npoints];
0d1b4cf8 2331
a23ba1c3 2332 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2333 Double_t lx=graph->GetX()[ipoint];
2334 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2335 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2336 if (index0<0) index0=0;
2337 if (index1>=npoints-1) index1=npoints-1;
2338 if ((index1-index0)>1){
2339 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2340 }else{
2341 outy[ipoint]=graph->GetY()[ipoint];
2342 }
2343 }
0d1b4cf8 2344 // TLinearFitter fitter(3,"pol2");
2345// for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2346// Double_t lx=graph->GetX()[ipoint];
2347// Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2348// Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2349// if (index0<0) index0=0;
2350// if (index1>=npoints-1) index1=npoints-1;
2351// fitter.ClearPoints();
2352// for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2353// if ((index1-index0)>1){
2354// outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2355// }else{
2356// outy[ipoint]=graph->GetY()[ipoint];
2357// }
2358// }
2359
2360
2361
a23ba1c3 2362 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2363 graph->GetY()[ipoint] = outy[ipoint];
2364 }
2365 delete[] outy;
2366}
2367
8166d5d7 2368Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
a23ba1c3 2369 //
2370 // Use constant interpolation outside of range
2371 //
2372 if (!graph) {
ef1807ef 2373 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
a23ba1c3 2374 return 0;
2375 }
5647625c 2376
a23ba1c3 2377 if (graph->GetN()<1){
ef1807ef 2378 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
a23ba1c3 2379 return 0;
2380 }
5647625c 2381
2382
a23ba1c3 2383 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2384 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
5647625c 2385
ef1807ef 2386 // AliInfo(Form("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref)));
5647625c 2387
2388 if(graph->GetN()==1)
2389 return graph->Eval(graph->GetX()[0]);
2390
2391
2392 return graph->Eval(xref);
a23ba1c3 2393}
2394
108a67db 2395Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2396 //
2397 // Use constant interpolation outside of range also for spline fits
2398 //
2399 if (!graph) {
ef1807ef 2400 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
108a67db 2401 return 0;
2402 }
2403 if (graph->GetKnots()<1){
ef1807ef 2404 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph");
108a67db 2405 return 0;
2406 }
2407 if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2408 if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2409 return graph->Eval( xref);
2410}
2411
1e722a63 2412Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2413 //
2414 // Filter DCS sensor information
2415 // ymin - minimal value
2416 // ymax - max value
2417 // maxdy - maximal deirivative
2418 // sigmaCut - cut on values and derivative in terms of RMS distribution
2419 // Return value - accepted fraction
2420 //
2421 // Algorithm:
2422 //
2423 // 0. Calculate median and rms of values in specified range
2424 // 1. Filter out outliers - median+-sigmaCut*rms
2425 // values replaced by median
2426 //
2427 AliSplineFit * fit = sensor->GetFit();
2428 if (!fit) return 0.;
2429 Int_t nknots = fit->GetKnots();
2430 if (nknots==0) {
2431 delete fit;
2432 sensor->SetFit(0);
2433 return 0;
2434 }
2435 //
2436 Double_t *yin0 = new Double_t[nknots];
2437 Double_t *yin1 = new Double_t[nknots];
2438 Int_t naccept=0;
2439
2440 for (Int_t iknot=0; iknot< nknots; iknot++){
2441 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2442 yin0[naccept] = fit->GetY0()[iknot];
2443 yin1[naccept] = fit->GetY1()[iknot];
2444 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2445 naccept++;
2446 }
2447 }
2448 if (naccept<1) {
2449 delete fit;
2450 sensor->SetFit(0);
a9f487a1 2451 delete [] yin0;
2452 delete [] yin1;
1e722a63 2453 return 0.;
2454 }
0d1b4cf8 2455
1e722a63 2456 Double_t medianY0=0, medianY1=0;
2457 Double_t rmsY0 =0, rmsY1=0;
2458 medianY0 = TMath::Median(naccept, yin0);
2459 medianY1 = TMath::Median(naccept, yin1);
2460 rmsY0 = TMath::RMS(naccept, yin0);
2461 rmsY1 = TMath::RMS(naccept, yin1);
2462 naccept=0;
2463 //
2464 // 1. Filter out outliers - median+-sigmaCut*rms
2465 // values replaced by median
2466 // if replaced the derivative set to 0
2467 //
2468 for (Int_t iknot=0; iknot< nknots; iknot++){
2469 Bool_t isOK=kTRUE;
2470 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2471 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2472 if (nknots<2) fit->GetY1()[iknot]=0;
2473 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2474 if (!isOK){
2475 fit->GetY0()[iknot]=medianY0;
2476 fit->GetY1()[iknot]=0;
2477 }else{
2478 naccept++;
2479 }
2480 }
2481 delete [] yin0;
2482 delete [] yin1;
2483 return Float_t(naccept)/Float_t(nknots);
2484}
2485
2486Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2487 //
2488 // Filter temperature array
2489 // tempArray - array of temperatures -
2490 // ymin - minimal accepted temperature - default 15
2491 // ymax - maximal accepted temperature - default 22
2492 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2493 // return value - fraction of filtered sensors
2494 const Double_t kMaxDy=0.1;
2495 Int_t nsensors=tempArray->NumSensors();
2496 if (nsensors==0) return 0.;
2497 Int_t naccept=0;
2498 for (Int_t isensor=0; isensor<nsensors; isensor++){
2499 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2500 if (!sensor) continue;
1e722a63 2501 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2502 if (sensor->GetFit()==0){
0d1b4cf8 2503 //delete sensor;
1e722a63 2504 tempArray->RemoveSensorNum(isensor);
2505 }else{
2506 naccept++;
2507 }
2508 }
2509 return Float_t(naccept)/Float_t(nsensors);
2510}
2511
a23ba1c3 2512
8166d5d7 2513void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2514 //
2515 // Filter CE data
1e722a63 2516 // Input parameters:
2517 // deltaT - smoothing window (in seconds)
2518 // cutAbs - max distance of the time info to the median (in time bins)
2519 // cutSigma - max distance (in the RMS)
2520 // pcstream - optional debug streamer to store original and filtered info
2521 // Hardwired parameters:
2522 // kMinPoints =10; // minimal number of points to define the CE
2523 // kMinSectors=12; // minimal number of sectors to define sideCE
2524 // Algorithm:
2525 // 0. Filter almost emty graphs (kMinPoints=10)
2526 // 1. calculate median and RMS per side
2527 // 2. Filter graphs - in respect with side medians
2528 // - cutAbs and cutDelta used
2529 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2530 // 4. Calculate mean for A side and C side
2531 //
2532 const Int_t kMinPoints =10; // minimal number of points to define the CE
2533 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2534 const Int_t kMinTime =400; // minimal arrival time of CE
2535 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
a23ba1c3 2536 Double_t medianY=0;
1e722a63 2537 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
a23ba1c3 2538 if (!cearray) return;
1e722a63 2539 Double_t tmin=-1;
2540 Double_t tmax=-1;
2541 //
2542 //
a23ba1c3 2543 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2cd3a697 2544 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
a23ba1c3 2545 if ( tempMapCE && cavernPressureCE){
1e722a63 2546 //
c5dfa2dc 2547 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2548 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2549 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2550 Bool_t isOK=kTRUE;
1e722a63 2551 if (isOK) {
2552 // recalculate P/T correction map for time of the CE
2553 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2554 driftCalib->SetName("driftPTCE");
2555 driftCalib->SetTitle("driftPTCE");
2556 cearray->AddLast(driftCalib);
2557 }
a23ba1c3 2558 }
1e722a63 2559 //
2560 // 0. Filter almost emty graphs
2561 //
a23ba1c3 2562
1e722a63 2563 for (Int_t i=0; i<72;i++){
a23ba1c3 2564 TGraph *graph= (TGraph*)arrT->At(i);
29797f44 2565 if (!graph) continue;
2566 graph->Sort();
a23ba1c3 2567 if (graph->GetN()<kMinPoints){
2568 arrT->AddAt(0,i);
2569 delete graph; // delete empty graph
2570 continue;
2571 }
1e722a63 2572 if (tmin<0) tmin = graph->GetX()[0];
2573 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2574 //
2575 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2576 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2577 }
2578 //
2579 // 1. calculate median and RMS per side
2580 //
2581 TArrayF arrA(100000), arrC(100000);
2582 Int_t nA=0, nC=0;
2583 Double_t medianA=0, medianC=0;
2584 Double_t rmsA=0, rmsC=0;
2585 for (Int_t isec=0; isec<72;isec++){
2586 TGraph *graph= (TGraph*)arrT->At(isec);
2587 if (!graph) continue;
2588 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2589 if (graph->GetY()[ipoint]<kMinTime) continue;
2590 if (nA>=arrA.fN) arrA.Set(nA*2);
2591 if (nC>=arrC.fN) arrC.Set(nC*2);
2592 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2593 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2594 }
2595 }
2596 if (nA>0){
2597 medianA=TMath::Median(nA,arrA.fArray);
2598 rmsA =TMath::RMS(nA,arrA.fArray);
2599 }
2600 if (nC>0){
2601 medianC=TMath::Median(nC,arrC.fArray);
2602 rmsC =TMath::RMS(nC,arrC.fArray);
2603 }
2604 //
2605 // 2. Filter graphs - in respect with side medians
2606 //
2607 TArrayD vecX(100000), vecY(100000);
2608 for (Int_t isec=0; isec<72;isec++){
2609 TGraph *graph= (TGraph*)arrT->At(isec);
2610 if (!graph) continue;
2611 Double_t median = (isec%36<18) ? medianA: medianC;
2612 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2613 Int_t naccept=0;
29797f44 2614 // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2615 for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
1e722a63 2616 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2617 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2618 vecX[naccept]= graph->GetX()[ipoint];
2619 vecY[naccept]= graph->GetY()[ipoint];
2620 naccept++;
2621 }
2622 if (naccept<kMinPoints){
2623 arrT->AddAt(0,isec);
2624 delete graph; // delete empty graph
2625 continue;
2626 }
2627 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2628 delete graph;
2629 arrT->AddAt(graph2,isec);
2630 }
2631 //
2632 // 3. Cut in respect wit the graph median
2633 //
2634 for (Int_t i=0; i<72;i++){
2635 TGraph *graph= (TGraph*)arrT->At(i);
2636 if (!graph) continue;
2637 //
2638 // filter in range
2639 //
2640 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
a23ba1c3 2641 if (!graphTS0) continue;
2642 if (graphTS0->GetN()<kMinPoints) {
2643 delete graphTS0;
2644 delete graph;
2645 arrT->AddAt(0,i);
2646 continue;
2647 }
1e722a63 2648 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
f03e487b 2649 if (!graphTS) continue;
a23ba1c3 2650 graphTS->Sort();
2651 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2652 if (pcstream){
2653 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2654 (*pcstream)<<"filterCE"<<
2655 "run="<<run<<
2656 "isec="<<i<<
2657 "mY="<<medianY<<
2658 "graph.="<<graph<<
2659 "graphTS0.="<<graphTS0<<
2660 "graphTS.="<<graphTS<<
2661 "\n";
2662 }
2663 delete graphTS0;
a23ba1c3 2664 arrT->AddAt(graphTS,i);
2665 delete graph;
2666 }
1e722a63 2667 //
2668 // Recalculate the mean time A side C side
2669 //
2670 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2671 Int_t meanPoints=(nA+nC)/72; // mean number of points
2672 for (Int_t itime=0; itime<200; itime++){
2673 nA=0, nC=0;
2674 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2675 for (Int_t i=0; i<72;i++){
2676 TGraph *graph= (TGraph*)arrT->At(i);
2677 if (!graph) continue;
2678 if (graph->GetN()<(meanPoints/4)) continue;
2679 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2680 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2681 }
2682 xA[itime]=time;
2683 xC[itime]=time;
2684 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2685 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2686 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2687 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2688 }
2689 //
2690 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2691 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2692 if (pcstream){
2693 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2694 (*pcstream)<<"filterAC"<<
2695 "run="<<run<<
2696 "nA="<<nA<<
2697 "nC="<<nC<<
2698 "rmsTA="<<rmsTA<<
2699 "rmsTC="<<rmsTC<<
2700 "\n";
2701 }
2702 //
2703 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2704 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2705 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2706 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2707 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2708 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2709 delete grA;
2710 delete grC;
2711 if (nA<kMinSectors) arrT->AddAt(0,72);
2712 else arrT->AddAt(graphTSA,72);
2713 if (nC<kMinSectors) arrT->AddAt(0,73);
2714 else arrT->AddAt(graphTSC,73);
a23ba1c3 2715}
2716
2717
8166d5d7 2718void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
a23ba1c3 2719 //
2720 // Filter Drift velocity measurement using the tracks
2721 // 0. remove outlyers - error based
2722 // cutSigma
2723 //
2724 //
2725 const Int_t kMinPoints=1; // minimal number of points to define value
2726 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2727 Double_t medianY=0;
2728 if (!arrT) return;
2729 for (Int_t i=0; i<arrT->GetEntries();i++){
14301155 2730 TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
a23ba1c3 2731 if (!graph) continue;
2732 if (graph->GetN()<kMinPoints){
2733 delete graph;
2734 arrT->AddAt(0,i);
2735 continue;
2736 }
5647625c 2737 TGraphErrors *graph2 = NULL;
2738 if(graph->GetN()<10) {
2739 graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2740 if (!graph2) {
2741 delete graph; arrT->AddAt(0,i); continue;
2742 }
2743 }
2744 else {
2745 graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2746 if (!graph2) {
2747 delete graph; arrT->AddAt(0,i); continue;
2748 }
a23ba1c3 2749 }
2750 if (graph2->GetN()<1) {
2751 delete graph; arrT->AddAt(0,i); continue;
2752 }
2753 graph2->SetName(graph->GetName());
2754 graph2->SetTitle(graph->GetTitle());
2755 arrT->AddAt(graph2,i);
2756 if (pcstream){
2757 (*pcstream)<<"filterTracks"<<
2758 "run="<<run<<
2759 "isec="<<i<<
2760 "mY="<<medianY<<
2761 "graph.="<<graph<<
2762 "graph2.="<<graph2<<
2763 "\n";
2764 }
2765 delete graph;
2766 }
2767}
2768
2769
a23ba1c3 2770
2771
2772
2773Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2774 //
2775 //
2776 // get laser time offset
2777 // median around timeStamp+-deltaT
2778 // QA - chi2 needed for later usage - to be added
2779 // - currently cut on error
2780 //
2781 Int_t kMinPoints=1;
2782 Double_t kMinDelay=0.01;
2783 Double_t kMinDelayErr=0.0001;
2784
2785 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2786 if (!array) return 0;
2787 TGraphErrors *tlaser=0;
2788 if (array){
2789 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2790 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2791 }
2792 if (!tlaser) return 0;
2793 Int_t npoints0= tlaser->GetN();
2794 if (npoints0==0) return 0;
2795 Double_t *xlaser = new Double_t[npoints0];
2796 Double_t *ylaser = new Double_t[npoints0];
2797 Int_t npoints=0;
2798 for (Int_t i=0;i<npoints0;i++){
2799 //printf("%d\n",i);
2800 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2801 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2802 xlaser[npoints]=tlaser->GetX()[npoints];
2803 ylaser[npoints]=tlaser->GetY()[npoints];
2804 npoints++;
2805 }
2806 //
2807 //
2808 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2809 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2810 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2811 if (index0<0) index0=0;
2812 if (index1>=npoints-1) index1=npoints-1;
f03e487b 2813 if (index1-index0<kMinPoints) {
2814 delete [] ylaser;
2815 delete [] xlaser;
2816 return 0;
2817 }
a23ba1c3 2818 //
2819 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2820 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2821 delete [] ylaser;
2822 delete [] xlaser;
2823 return mean;
2824}
0d1b4cf8 2825
2826
2827
2828
8166d5d7 2829void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
0d1b4cf8 2830 //
2831 // Filter Goofie data
1fabc823 2832 // goofieArray - points will be filtered
2833 // deltaT - smmothing time window
2834 // cutSigma - outler sigma cut in rms
2835 // minVn, maxVd- range absolute cut for variable vd/pt
2836 // - to be tuned
0d1b4cf8 2837 //
0d1b4cf8 2838 // Ignore goofie if not enough points
2839 //
2840 const Int_t kMinPoints = 3;
2841 //
2842
2843 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2844 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2845 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2846 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2847 if (!graphvd) return;
2848 if (graphvd->GetN()<kMinPoints){
2849 delete graphvd;
2850 goofieArray->GetSensorNum(2)->SetGraph(0);
2851 return;
2852 }
2853 //
2854 // 1. Caluclate medians of critical variables
2855 // drift velcocity
2856 // P/T
2857 // area near peak
2858 // area far peak
2859 //
2860 Double_t medianpt=0;
1fabc823 2861 Double_t medianvd=0, sigmavd=0;
0d1b4cf8 2862 Double_t medianan=0;
2863 Double_t medianaf=0;
1fabc823 2864 Int_t entries=graphvd->GetN();
2865 Double_t yvdn[10000];
2866 Int_t nvd=0;
2867 //
2868 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2869 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2870 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2871 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2872 yvdn[nvd++]=graphvd->GetY()[ipoint];
2873 }
2874 if (nvd<kMinPoints){
2875 delete graphvd;
2876 goofieArray->GetSensorNum(2)->SetGraph(0);
2877 return;
2878 }
2879 //
2880 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2881 if (nuni>=kMinPoints){
2882 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2883 }else{
2884 medianvd = TMath::Median(nvd, yvdn);
2885 }
2886
0d1b4cf8 2887 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2888 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2889 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2890 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2891 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2892 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
0d1b4cf8 2893 delete graphpt0;
2894 delete graphpt1;
2895 delete graphan0;
2896 delete graphan1;
2897 delete graphaf0;
2898 delete graphaf1;
2899 //
2900 // 2. Make outlyer graph
2901 //
1fabc823 2902 Int_t nOK=0;
0d1b4cf8 2903 TGraph graphOut(*graphvd);
2904 for (Int_t i=0; i<entries;i++){
2905 //
2906 Bool_t isOut=kFALSE;
1fabc823 2907 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2908 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2909
2910 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2911 isOut|=kTRUE;
0d1b4cf8 2912 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
1fabc823 2913 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2914 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
0d1b4cf8 2915 graphOut.GetY()[i]= (isOut)?1:0;
1fabc823 2916 if (!isOut) nOK++;
0d1b4cf8 2917 }
1fabc823 2918 if (nOK<kMinPoints) {
2919 delete graphvd;
2920 goofieArray->GetSensorNum(2)->SetGraph(0);
2921 return;
2922 }
0d1b4cf8 2923 //
2924 // 3. Filter out outlyers - and smooth
2925 //
2926 TVectorF vmedianArray(goofieArray->NumSensors());
2927 TVectorF vrmsArray(goofieArray->NumSensors());
2928 Double_t xnew[10000];
2929 Double_t ynew[10000];
2930 TObjArray junk;
2931 junk.SetOwner(kTRUE);
2932 Bool_t isOK=kTRUE;
2933 //
2934 //
2935 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2936 isOK=kTRUE;
2937 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2938 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2939 //
2940 if (!sensor) continue;
2941 graphOld = sensor->GetGraph();
2942 if (graphOld) {
2943 sensor->SetGraph(0);
2944 Int_t nused=0;
2945 for (Int_t i=0;i<entries;i++){
2946 if (graphOut.GetY()[i]>0.5) continue;
2947 xnew[nused]=graphOld->GetX()[i];
2948 ynew[nused]=graphOld->GetY()[i];
2949 nused++;
2950 }
2951 graphNew = new TGraph(nused,xnew,ynew);
2952 junk.AddLast(graphNew);
2953 junk.AddLast(graphOld);
2954 Double_t median=0;
2955 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2956 if (graphNew0!=0){
2957 junk.AddLast(graphNew0);
2958 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2959 if (graphNew1!=0){
1fabc823 2960 junk.AddLast(graphNew1);
0d1b4cf8 2961 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2962 if (graphNew2!=0) {
1fabc823 2963 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
0d1b4cf8 2964 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2965 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2966 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
ef1807ef 2967 // AliInfo(Form("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]));
0d1b4cf8 2968 vmedianArray[isensor]=median;
0d1b4cf8 2969 //
2970 }
2971 }
2972 }
2973 }
2974 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2975 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2976 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2977 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2978 (*pcstream)<<"goofieA"<<
2979 Form("isOK_%d.=",isensor)<<isOK<<
2980 Form("s_%d.=",isensor)<<sensor<<
2981 Form("gr_%d.=",isensor)<<graphOld<<
2982 Form("gr0_%d.=",isensor)<<graphNew0<<
2983 Form("gr1_%d.=",isensor)<<graphNew1<<
2984 Form("gr2_%d.=",isensor)<<graphNew2;
1fabc823 2985 if (isOK) sensor->SetGraph(graphNew2);
2986 }
2987 (*pcstream)<<"goofieA"<<
2988 "vmed.="<<&vmedianArray<<
2989 "vrms.="<<&vrmsArray<<
2990 "\n";
2991 junk.Delete(); // delete temoprary graphs
0d1b4cf8 2992
2993}
2994
2995
a980538f 2996
2997
2998
8166d5d7 2999TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
a980538f 3000 //
3001 // Make a statistic matrix
3002 // Input parameters:
3003 // array - TObjArray of AliRelKalmanAlign
3004 // minFraction - minimal ration of accepted tracks
3005 // minStat - minimal statistic (number of accepted tracks)
3006 // maxvd - maximal deviation for the 1
3007 // Output matrix:
3008 // columns - Mean, Median, RMS
3009 // row - parameter type (rotation[3], translation[3], drift[3])
3010 if (!array) return 0;
3011 if (array->GetEntries()<=0) return 0;
cc65e4f5 3012 // Int_t entries = array->GetEntries();
a980538f 3013 Int_t entriesFast = array->GetEntriesFast();
3014 TVectorD state(9);
3015 TVectorD *valArray[9];
3016 for (Int_t i=0; i<9; i++){
3017 valArray[i] = new TVectorD(entriesFast);
3018 }
3019 Int_t naccept=0;
3020 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3021 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3022 if (!kalman) continue;
3023 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3024 if (kalman->GetNUpdates()<minStat) continue;
938ce618 3025 if (Float_t(kalman->GetNUpdates())/Float_t(kalman->GetNTracks())<minFraction) continue;
a980538f 3026 kalman->GetState(state);
3027 for (Int_t ipar=0; ipar<9; ipar++)
3028 (*valArray[ipar])[naccept]=state[ipar];
3029 naccept++;
3030 }
5647625c 3031 //if (naccept<2) return 0;
3032 if (naccept<1) return 0;
a980538f 3033 TMatrixD *pstat=new TMatrixD(9,3);
3034 TMatrixD &stat=*pstat;
3035 for (Int_t ipar=0; ipar<9; ipar++){
3036 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3037 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3038 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3039 }
177f532f 3040 for (Int_t i=0; i<9; i++){
3041 delete valArray[i];
3042 }
a980538f 3043 return pstat;
3044}
3045
3046
8166d5d7 3047TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
a980538f 3048 //
3049 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3050 // Input:
3051 // array - input array
3052 // stat - mean parameters statistic
3053 // direction -
3054 // sigmaCut - maximal allowed deviation from mean in terms of RMS
3055 if (!array) return 0;
3056 if (array->GetEntries()<=0) return 0;
32100b1d 3057 if (!(&stat)) return 0;
cc65e4f5 3058 // error increase in 1 hour
3059 const Double_t kerrsTime[9]={
3060 0.00001, 0.00001, 0.00001,
3061 0.001, 0.001, 0.001,
dfb83639 3062 0.002, 0.01, 0.001};
cc65e4f5 3063 //
3064 //
a980538f 3065 Int_t entries = array->GetEntriesFast();
3066 TObjArray *sArray= new TObjArray(entries);
3067 AliRelAlignerKalman * sKalman =0;
3068 TVectorD state(9);
3069 for (Int_t i=0; i<entries; i++){
3070 Int_t index=(direction)? entries-i-1:i;
3071 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3072 if (!kalman) continue;
3073 Bool_t isOK=kTRUE;
3074 kalman->GetState(state);
3075 for (Int_t ipar=0; ipar<9; ipar++){
3076 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3077 }
3078 if (!sKalman &&isOK) {
3079 sKalman=new AliRelAlignerKalman(*kalman);
3080 sKalman->SetRejectOutliers(kFALSE);
3081 sKalman->SetRunNumber(kalman->GetRunNumber());
3082 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3083 }
3084 if (!sKalman) continue;
3085 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
cc65e4f5 3086 for (Int_t ipar=0; ipar<9; ipar++){
3087// (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3088// (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3089// (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3090 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3091 }
a980538f 3092 sKalman->SetRunNumber(kalman->GetRunNumber());
3093 if (!isOK) sKalman->SetRunNumber(0);
3094 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3095 if (!isOK) continue;
3096 sKalman->SetRejectOutliers(kFALSE);
3097 sKalman->SetRunNumber(kalman->GetRunNumber());
3098 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3099 sKalman->Merge(kalman);
3100 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3101 //sKalman->Print();
3102 }
3103 return sArray;
3104}
3105
8166d5d7 3106TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
cc65e4f5 3107 //
3108 // Merge 2 RelKalman arrays
3109 // Input:
3110 // arrayP - rel kalman in direction plus
3111 // arrayM - rel kalman in direction minus
3112 if (!arrayP) return 0;
3113 if (arrayP->GetEntries()<=0) return 0;
3114 if (!arrayM) return 0;
3115 if (arrayM->GetEntries()<=0) return 0;
3a4b0aed 3116
cc65e4f5 3117 Int_t entries = arrayP->GetEntriesFast();
3118 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3a4b0aed 3119
cc65e4f5 3120 for (Int_t i=0; i<entries; i++){
3121 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3122 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3123 if (!kalmanP) continue;
3124 if (!kalmanM) continue;
3a4b0aed 3125
3126 AliRelAlignerKalman *kalman = NULL;
3127 if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3128 kalman = new AliRelAlignerKalman(*kalmanP);
3129 kalman->Merge(kalmanM);
3130 }
3131 else if (kalmanP->GetRunNumber() == 0) {
3132 kalman = new AliRelAlignerKalman(*kalmanM);
3133 }
3134 else if (kalmanM->GetRunNumber() == 0) {
3135 kalman = new AliRelAlignerKalman(*kalmanP);
3136 }
3137 else
3138 continue;
3139
cc65e4f5 3140 array->AddAt(kalman,i);
3141 }
3a4b0aed 3142
cc65e4f5 3143 return array;
3144}
c4e4498c 3145
3146//_____________________________________________________________________________________
3147TTree* AliTPCcalibDButil::ConnectGainTrees(TString baseDir)
3148{
3149 //
3150 // baseDir: Base directory with the raw Kr calibration trees
3151 // and the trees from the calibQA
3152 // it assumes to following structure below:
3153 // KryptonCalib/<year>/calibKr/calibKr.<year>.<id>.root
3154 // calibQAdEdx/<year>/calibQA.<year>.<perid>.tree.root
3155 // map/treeMapping.root
3156 //
3157
3158
3159 // === add main tree, which will be a mapping file ================
3160 TFile *fin = TFile::Open(Form("%s/map/treeMapping.root",baseDir.Data()));
3161 gROOT->cd();
3162 TTree *tMain = (TTree*)fin->Get("calPads");
a08e21e8 3163 tMain->SetAlias("Pad0","MapPadLength.fElements==7.5"); // pad types
3164 tMain->SetAlias("Pad1","MapPadLength.fElements==10.0");
3165 tMain->SetAlias("Pad2","MapPadLength.fElements==15.0");
3166 //
3167 tMain->SetAlias("dRowNorm0","(row.fElements/32-1)"); // normalized distance to the center of the chamber in lx (-1,1)
3168 tMain->SetAlias("dRowNorm1","(row.fElements/32-1)"); //
3169 tMain->SetAlias("dRowNorm2","((row.fElements-63)/16-1)"); //
3170 tMain->SetAlias("dRowNorm","(row.fElements<63)*(row.fElements/32-1)+(row.fElements>63)*((row.fElements-63)/16-1)");
3171 tMain->SetAlias("dPhiNorm","(ly.fElements/lx.fElements)/(pi/18)"); // normalized distance to the center of the chamber in phi (-1,1)
3172 //
3173 tMain->SetAlias("fsector","(sector%36)+(0.5+9*atan2(ly.fElements,lx.fElements)/pi)"); // float sector number
3174 tMain->SetAlias("posEdge","((pi/18.)-abs(atan(ly.fElements/lx.fElements)))*lx.fElements"); //distance to the edge
c4e4498c 3175
3176 // === add the krypton calibration trees ==========================
3177 TString inputTreesKrCalib = gSystem->GetFromPipe(Form("ls %s/KryptonCalib/20*/calibKr/*.tree.root",baseDir.Data()));
3178 TObjArray *arrInputTreesKrCalib = inputTreesKrCalib.Tokenize("\n");
3179
3180 for (Int_t itree=0; itree<arrInputTreesKrCalib->GetEntriesFast(); ++itree) {
3181 TFile *fin2 = TFile::Open(arrInputTreesKrCalib->At(itree)->GetName());
3182 TTree *tin = (TTree*)fin2->Get("calPads");
3183 gROOT->cd();
3184 TString friendName=gSystem->BaseName(arrInputTreesKrCalib->At(itree)->GetName());
3185 friendName.ReplaceAll("calibKr.","");
3186 friendName.ReplaceAll(".tree.root","");
a08e21e8 3187 friendName="Kr."+friendName;
c4e4498c 3188 tMain->AddFriend(tin,friendName.Data());
3189
3190 // set aliases
3191
3192 // TODO: finish implementation of alias via lists
3193// const Int_t nbranchAlias = 2;
3194// const char* branchNames[nbranchAlias]={"spectrMean","fitMean"};
3195// const Int_t nbranchStat = 2;
3196// const char* statNames[nbranchStat] = {"Median","LTM"};
3197//
3198// for (Int_t iname=0; iname<nbranchAlias; ++iname) {
3199// TString branchName = TString::Format("%s.%s", friendName.Data(), bNames[i]);
3200//
3201// for (Int_t istat=0; istat<nbranchStat; ++istat) {
3202//
3203// }
3204// }
3205
3206 tMain->SetAlias((friendName+".spectrMean_LTMRatio").Data(),
05ff31fd 3207 TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_LTM+0))",
c4e4498c 3208 friendName.Data(),friendName.Data()).Data());
3209
3210 tMain->SetAlias((friendName+".spectrMean_MedianRatio").Data(),
05ff31fd 3211 TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_Median+0))",
c4e4498c 3212 friendName.Data(),friendName.Data()).Data());
a08e21e8 3213 tMain->SetAlias((friendName+".spectrMean_MeanRatio").Data(),
05ff31fd 3214 TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_Mean+0))",
a08e21e8 3215 friendName.Data(),friendName.Data()).Data());
e2179c10 3216 tMain->SetAlias((friendName+".spectrMean_MeanRatio").Data(),
3217 TString::Format("(%s.spectrMean.fElements/%s.spectrMean_Mean)",
3218 friendName.Data(),friendName.Data()).Data());
c4e4498c 3219
3220 tMain->SetAlias((friendName+".fitMean_LTMRatio").Data(),
05ff31fd 3221 TString::Format("(%s.fitMean.fElements/(%s.fitMean_LTM+0))",
3222 friendName.Data(),friendName.Data()).Data());
3223 tMain->SetAlias((friendName+".fitRMS_LTMRatio").Data(),
3224 TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_LTM+0))",
c4e4498c 3225 friendName.Data(),friendName.Data()).Data());
3226
3227 tMain->SetAlias((friendName+".fitMean_MedianRatio").Data(),
05ff31fd 3228 TString::Format("(%s.fitMean.fElements/(%s.fitMean_Median+0))",
3229 friendName.Data(),friendName.Data()).Data());
3230 tMain->SetAlias((friendName+".fitRMS_MedianRatio").Data(),
3231 TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_Median+0))",
c4e4498c 3232 friendName.Data(),friendName.Data()).Data());
a08e21e8 3233 tMain->SetAlias((friendName+".fitMean_MeanRatio").Data(),
05ff31fd 3234 TString::Format("(%s.fitMean.fElements/(%s.fitMean_Mean+0))",
3235 friendName.Data(),friendName.Data()).Data());
3236 tMain->SetAlias((friendName+".fitRMS_MeanRatio").Data(),
3237 TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_Mean+0))",
a08e21e8 3238 friendName.Data(),friendName.Data()).Data());
e2179c10 3239 tMain->SetAlias((friendName+".fitMean_MeanRatio").Data(),
3240 TString::Format("(%s.fitMean.fElements/%s.fitMean_Mean)",
3241 friendName.Data(),friendName.Data()).Data());
c4e4498c 3242
3243 }
3244
3245 // === add the calibQA trees ======================================
3246 TString inputTreesQACalib = gSystem->GetFromPipe(Form("ls %s/calibQAdEdx/20*/*.tree.root",baseDir.Data()));
3247 TObjArray *arrInputTreesQACalib = inputTreesQACalib.Tokenize("\n");
3248
3249 for (Int_t itree=0; itree<arrInputTreesQACalib->GetEntriesFast(); ++itree) {
3250 TFile *fin2 = TFile::Open(arrInputTreesQACalib->At(itree)->GetName());
3251 TTree *tin = (TTree*)fin2->Get("calPads");
3252 gROOT->cd();
3253 TString friendName=gSystem->BaseName(arrInputTreesQACalib->At(itree)->GetName());
3254 friendName.ReplaceAll("calibQA.","");
3255 friendName.ReplaceAll(".tree.root","");
a08e21e8 3256 friendName="QA."+friendName;
c4e4498c 3257 tMain->AddFriend(tin,friendName.Data());
3258
3259 // set aliases
3260 tMain->SetAlias((friendName+".MaxCharge_LTMRatio").Data(),
3261 TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_LTM)",
3262 friendName.Data(),friendName.Data()).Data());
3263
3264 tMain->SetAlias((friendName+".MaxCharge_MedianRatio").Data(),
3265 TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_Median)",
3266 friendName.Data(),friendName.Data()).Data());
a08e21e8 3267 tMain->SetAlias((friendName+".MaxCharge_MeanRatio").Data(),
3268 TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_Mean)",
3269 friendName.Data(),friendName.Data()).Data());
c4e4498c 3270
3271 tMain->SetAlias((friendName+".MeanCharge_LTMRatio").Data(),
3272 TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_LTM)",
3273 friendName.Data(),friendName.Data()).Data());
3274
3275 tMain->SetAlias((friendName+".MeanCharge_MedianRatio").Data(),
3276 TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_Median)",
3277 friendName.Data(),friendName.Data()).Data());
a08e21e8 3278 tMain->SetAlias((friendName+".MeanCharge_MeanRatio").Data(),
3279 TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_Mean)",
3280 friendName.Data(),friendName.Data()).Data());
c4e4498c 3281
3282 }
3283
3284 return tMain;
3285}
3286
06a262f6 3287
3288//_____________________________________________________________________________________
3289TTree* AliTPCcalibDButil::ConnectPulserTrees(TString baseDir, TTree *tMain)
3290{
3291 //
3292 // baseDir: Base directory with Pulser information
3293 // TTrees are added to the base tree as a friend tree
3294 //
3295 // === add the calibPulser trees ======================================
3296 TString inputTreesPulserCalib = gSystem->GetFromPipe(Form("ls %s/calibPulser/20*/*.tree.root",baseDir.Data()));
3297 TObjArray *arrInputTreesPulserCalib = inputTreesPulserCalib.Tokenize("\n");
3298 for (Int_t itree=0; itree<arrInputTreesPulserCalib->GetEntriesFast(); ++itree) {
3299 TFile *fin2 = TFile::Open(arrInputTreesPulserCalib->At(itree)->GetName());
3300 TTree *tin = (TTree*)fin2->Get("calPads");
3301 gROOT->cd();
3302 TString friendName=gSystem->BaseName(arrInputTreesPulserCalib->At(itree)->GetName());
3303 friendName.ReplaceAll("calibPulser.","");
3304 friendName.ReplaceAll(".tree.root","");
3305 friendName="Pulser."+friendName;
3306 tMain->AddFriend(tin,friendName.Data());
3307 // set aliases
3308
3309 tMain->SetAlias((friendName+".CEQmean_LTMRatio").Data(),
3310 TString::Format("(%s.CEQmean.fElements/%s.CEQmean_LTM)",
3311 friendName.Data(),friendName.Data()).Data());
3312 tMain->SetAlias((friendName+".CEQmean_MedianRatio").Data(),
3313 TString::Format("(%s.CEQmean.fElements/%s.CEQmean_Median)",
3314 friendName.Data(),friendName.Data()).Data());
3315 tMain->SetAlias((friendName+".CEQmean_MeanRatio").Data(),
3316 TString::Format("(%s.CEQmean.fElements/%s.CEQmean_Mean)",
3317 friendName.Data(),friendName.Data()).Data());
3318 //
3319 tMain->SetAlias((friendName+".CETmean_LTMDelta").Data(),
3320 TString::Format("(%s.CETmean.fElements-%s.CETmean_LTM)",
3321 friendName.Data(),friendName.Data()).Data());
3322 tMain->SetAlias((friendName+".CETmean_MedianDelta").Data(),
c4f17a0d 3323 TString::Format("(%s.CETmean.fElements-%s.CETmean_Median)",
06a262f6 3324 friendName.Data(),friendName.Data()).Data());
3325 tMain->SetAlias((friendName+".CETmean_MeanDelta").Data(),
3326 TString::Format("(%s.CETmean.fElements-%s.CETmean_Mean)",
3327 friendName.Data(),friendName.Data()).Data());
3328 }
c4f17a0d 3329 return tMain;
06a262f6 3330}
3331
c4f17a0d 3332
525f22d2 3333TTree* AliTPCcalibDButil::ConnectDistortionTrees(TString baseDir, TString selection, TTree *tMain){
c4f17a0d 3334 //
3335 // baseDir: Base directory with Distortion information
3336 // TTrees are added to the base tree as a friend tree
3337 // If base tree not provide - first tree from list is used as base
3338 //
3339 // === add the calibDistortion trees ======================================
3340 //TString inputTreesDistortionCalib = gSystem->GetFromPipe(Form("ls %s/calibDistortion/20*/*.tree.root",baseDir.Data()));
3341 // TString baseDir="$NOTES/reconstruction/distortionFit/"; TTree *tMain=0;
525f22d2 3342 // AliTPCcalibDButil::ConnectDistortionTrees("$NOTES/reconstruction/distortionFit/", "calibTimeResHisto.root", 0);
c4f17a0d 3343
525f22d2 3344 TString inputTreesDistortionCalib = "";
3345 if (selection.Contains(".list")){
3346 inputTreesDistortionCalib=gSystem->GetFromPipe(Form("cat %s",selection.Data()));
3347 }else{
3348 inputTreesDistortionCalib=gSystem->GetFromPipe(Form("find %s -iname \"%s\"",baseDir.Data(),selection.Data()));
3349 }
c4f17a0d 3350 TObjArray *arrInputTreesDistortionCalib = inputTreesDistortionCalib.Tokenize("\n");
3351 //
3352 for (Int_t itree=0; itree<arrInputTreesDistortionCalib->GetEntriesFast(); ++itree) {
3353 TFile *finput= TFile::Open(arrInputTreesDistortionCalib->At(itree)->GetName());
3354 TString strFile=arrInputTreesDistortionCalib->At(itree)->GetName();
3355 TObjArray *path=strFile.Tokenize("/");
3356 Int_t plength=path->GetEntries();
3357 if (!finput) continue;
3358 TList* list = finput->GetListOfKeys();
3359 Int_t nkeys=list->GetEntries();
3360 for (Int_t ikey=0; ikey<nkeys; ikey++){
3361 TKey * key = (TKey*)list->At(ikey);
3362 if (strstr(key->GetClassName(),"TTree")==0) continue;
3363 TTree * tree = dynamic_cast<TTree*>(finput->Get(list->At(ikey)->GetName()));
3364 if (!tree) continue;
3365 TString friendName=TString::Format("%s.%s.%s",path->At(plength-3)->GetName(),path->At(plength-2)->GetName(), tree->GetName());
3366 ::Info("AliTPCcalibDButil::ConnectDistortionTrees","%s",friendName.Data());
3367 if (tMain==0) tMain=tree;
3368 tMain->AddFriend(tree,friendName.Data());
3369 }
3370 }
3371 // tMain->SetAlias("");
3372 return tMain;
3373}
3374
3375
3376//_____________________________________________________________________________________
085597e9 3377TTree* AliTPCcalibDButil::ConnectCalPadTrees(TString baseDir, TString pattern, TTree *tMain, Bool_t checkAliases)
c4f17a0d 3378{
3379 //
3380 // baseDir: Base directory with per Pad information
3381 // TTrees are added to the base tree as a friend tree
085597e9 3382 // Example usage
3383 // TString baseDir="/hera/alice/fsozzi/summarymaps/calib2/"; // prefix directory with calibration with slash at the end
3384 // TString pattern="QA/*/*root";
3385 // TTree * tree = AliTPCcalibDButil::ConnectCalPadTrees(baseDir,pattern,0); //create tree and attach calibration as friends
3386
c4f17a0d 3387 //
3388 // === add the calibPulser trees ======================================
3389 TString inputTreesCalPad = gSystem->GetFromPipe(Form("ls %s/%s",baseDir.Data(), pattern.Data()));
3390 TObjArray *arrInputTreesCalPad = inputTreesCalPad.Tokenize("\n");
3391 //
3392 for (Int_t itree=0; itree<arrInputTreesCalPad->GetEntriesFast(); ++itree) {
3393 TFile *fin2 = TFile::Open(arrInputTreesCalPad->At(itree)->GetName());
3394 TTree *tin = (TTree*)fin2->Get("calPads");
3395 gROOT->cd();
3396 TString friendName=arrInputTreesCalPad->At(itree)->GetName();
085597e9 3397 friendName.ReplaceAll("//","/");
c4f17a0d 3398 friendName.ReplaceAll(baseDir.Data(),"");
085597e9 3399 friendName.ReplaceAll("^/","");
c4f17a0d 3400 friendName.ReplaceAll("/",".");
3401 friendName.ReplaceAll(".root","");
c4f17a0d 3402 printf("%s\n", friendName.Data());
3403 ::Info("AliTPCcalibDButil::ConnectCalPadTrees","%s",friendName.Data());
3404 if (tMain==0) tMain=tin;
3405 tMain->AddFriend(tin,friendName.Data());
085597e9 3406 TObjArray * branches=tin->GetListOfBranches();
3407 Int_t nBranches=branches->GetEntries();
3408 for (Int_t ibranch=0; ibranch<nBranches; ibranch++){
3409 TString bname=branches->At(ibranch)->GetName();
3410 if (bname.Contains(".")>0){
3411 bname.ReplaceAll(".","");
3412 // replace elements
3413 tin->SetAlias((bname).Data(), (bname+".fElements").Data());
3414 tMain->SetAlias((friendName+"."+bname).Data(), (friendName+"."+bname+".fElements").Data());
3415 //
3416 // make normalized values per chamber
3417 //
3418 if (branches->FindObject(bname+"_LTM")!=0){
3419 tMain->SetAlias((friendName+"."+bname+"_MeanRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_Mean").Data());
3420 tMain->SetAlias((friendName+"."+bname+"_MedianRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_Median").Data());
3421 tMain->SetAlias((friendName+"."+bname+"_LTMRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_LTM").Data());
3422 tMain->SetAlias((friendName+"."+bname+"_MeanDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_Mean").Data());
3423 tMain->SetAlias((friendName+"."+bname+"_MedianDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_Median").Data());
3424 tMain->SetAlias((friendName+"."+bname+"_LTMDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_LTM").Data());
3425 }
8bc3ddf0 3426 // if (branches->FindObject(bname+"_Med3.")!=0){
3427 tMain->SetAlias((friendName+"."+bname+"_Med3Ratio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Med3.fElements").Data());
3428 tMain->SetAlias((friendName+"."+bname+"_Med5Ratio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Med5.fElements").Data());
3429 tMain->SetAlias((friendName+"."+bname+"_Par6GRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Par6G.fElements").Data());
3430 tMain->SetAlias((friendName+"."+bname+"_Med3Delta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Med3.fElements").Data());
3431 tMain->SetAlias((friendName+"."+bname+"_Med5Delta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Med5.fElements").Data());
3432 tMain->SetAlias((friendName+"."+bname+"_Par6GDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Par6G.fElements").Data());
3433 // }
085597e9 3434 }
3435 }
3436 }
3437 //
3438 //
3439 //
3440 if (checkAliases){
3441 // to be implemented
c4f17a0d 3442 }
3443 return tMain;
3444}
3445
3446
3447