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