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