TPC/AliTPCCalibCE.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
  20 //-------------------------------------------------------
  21 //          Implementation of the TPC pulser calibration
  22 //
  23 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
  24 //
  25 //
  26 //-------------------------------------------------------
  27
  28
  29 /* $Id$ */
  30
  31
  32
  33 //Root includes
  34 #include <TObjArray.h>
  35 #include <TH1F.h>
  36 #include <TH2S.h>
  37 #include <TString.h>
  38 #include <TVectorF.h>
  39 #include <TVectorD.h>
  40 #include <TMatrixD.h>
  41 #include <TMath.h>
  42 #include <TGraph.h>
  43
  44 #include <TDirectory.h>
  45 #include <TSystem.h>
  46 #include <TFile.h>
  47
  48 //AliRoot includes
  49 #include "AliRawReader.h"
  50 #include "AliRawReaderRoot.h"
  51 #include "AliRawReaderDate.h"
  52 #include "AliRawEventHeaderBase.h"
  53 #include "AliTPCRawStream.h"
  54 #include "AliTPCcalibDB.h"
  55 #include "AliTPCCalROC.h"
  56 #include "AliTPCCalPad.h"
  57 #include "AliTPCROC.h"
  58 #include "AliTPCParam.h"
  59 #include "AliTPCCalibCE.h"
  60 #include "AliMathBase.h"
  61 #include "TTreeStream.h"
  62
  63 //date
  64 #include "event.h"
  65 ClassImp(AliTPCCalibCE) /*FOLD00*/
  66
  67 AliTPCCalibCE::AliTPCCalibCE() : /*FOLD00*/
  68     TObject(),
  69     fFirstTimeBin(650),
  70     fLastTimeBin(1000),
  71     fNbinsT0(200),
  72     fXminT0(-5),
  73     fXmaxT0(5),
  74     fNbinsQ(200),
  75     fXminQ(1),
  76     fXmaxQ(40),
  77     fNbinsRMS(100),
  78     fXminRMS(0.1),
  79     fXmaxRMS(5.1),
  80     fLastSector(-1),
  81     fOldRCUformat(kTRUE),
  82     fROC(AliTPCROC::Instance()),
  83     fParam(new AliTPCParam),
  84     fPedestalTPC(0x0),
  85     fPadNoiseTPC(0x0),
  86     fPedestalROC(0x0),
  87     fPadNoiseROC(0x0),
  88     fCalRocArrayT0(72),
  89     fCalRocArrayQ(72),
  90     fCalRocArrayRMS(72),
  91     fCalRocArrayOutliers(72),
  92     fHistoQArray(72),
  93     fHistoT0Array(72),
  94     fHistoRMSArray(72),
  95     fHistoTmean(72),
  96     fParamArrayEvent(1000),
  97     fParamArrayEventPol1(72),
  98     fParamArrayEventPol2(72),
  99     fTMeanArrayEvent(1000),
 100     fQMeanArrayEvent(1000),
 101     fVEventTime(1000),
 102     fVEventNumber(1000),
 103     fNevents(0),
 104     fTimeStamp(0),
 105     fRunNumber(-1),
 106     fOldRunNumber(-1),
 107     fPadTimesArrayEvent(72),
 108     fPadQArrayEvent(72),
 109     fPadRMSArrayEvent(72),
 110     fPadPedestalArrayEvent(72),
 111     fCurrentChannel(-1),
 112     fCurrentSector(-1),
 113     fCurrentRow(-1),
 114     fMaxPadSignal(-1),
 115     fMaxTimeBin(-1),
 116     fPadSignal(1024),
 117     fPadPedestal(0),
 118     fPadNoise(0),
 119     fVTime0Offset(72),
 120     fVTime0OffsetCounter(72),
 121     fVMeanQ(72),
 122     fVMeanQCounter(72),
 123 //    fHTime0(0x0),
 124     fEvent(-1),
 125     fDebugStreamer(0x0),
 126     fDebugLevel(0)
 127 {
 128     //
 129     // AliTPCSignal default constructor
 130     //
 131 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
 132 }
 133 //_____________________________________________________________________
 134 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
 135     TObject(sig),
 136     fFirstTimeBin(sig.fFirstTimeBin),
 137     fLastTimeBin(sig.fLastTimeBin),
 138     fNbinsT0(sig.fNbinsT0),
 139     fXminT0(sig.fXminT0),
 140     fXmaxT0(sig.fXmaxT0),
 141     fNbinsQ(sig.fNbinsQ),
 142     fXminQ(sig.fXminQ),
 143     fXmaxQ(sig.fXmaxQ),
 144     fNbinsRMS(sig.fNbinsRMS),
 145     fXminRMS(sig.fXminRMS),
 146     fXmaxRMS(sig.fXmaxRMS),
 147     fLastSector(-1),
 148     fOldRCUformat(kTRUE),
 149     fROC(AliTPCROC::Instance()),
 150     fParam(new AliTPCParam),
 151     fPedestalTPC(0x0),
 152     fPadNoiseTPC(0x0),
 153     fPedestalROC(0x0),
 154     fPadNoiseROC(0x0),
 155     fCalRocArrayT0(72),
 156     fCalRocArrayQ(72),
 157     fCalRocArrayRMS(72),
 158     fCalRocArrayOutliers(72),
 159     fHistoQArray(72),
 160     fHistoT0Array(72),
 161     fHistoRMSArray(72),
 162     fHistoTmean(72),
 163     fParamArrayEvent(1000),
 164     fParamArrayEventPol1(72),
 165     fParamArrayEventPol2(72),
 166     fTMeanArrayEvent(1000),
 167     fQMeanArrayEvent(1000),
 168     fVEventTime(1000),
 169     fVEventNumber(1000),
 170     fNevents(sig.fNevents),
 171     fTimeStamp(0),
 172     fRunNumber(-1),
 173     fOldRunNumber(-1),
 174     fPadTimesArrayEvent(72),
 175     fPadQArrayEvent(72),
 176     fPadRMSArrayEvent(72),
 177     fPadPedestalArrayEvent(72),
 178     fCurrentChannel(-1),
 179     fCurrentSector(-1),
 180     fCurrentRow(-1),
 181     fMaxPadSignal(-1),
 182     fMaxTimeBin(-1),
 183     fPadSignal(1024),
 184     fPadPedestal(0),
 185     fPadNoise(0),
 186     fVTime0Offset(72),
 187     fVTime0OffsetCounter(72),
 188     fVMeanQ(72),
 189     fVMeanQCounter(72),
 190 //    fHTime0(0x0),
 191     fEvent(-1),
 192     fDebugStreamer(0x0),
 193     fDebugLevel(sig.fDebugLevel)
 194 {
 195     //
 196     // AliTPCSignal default constructor
 197     //
 198
 199     for (Int_t iSec = 0; iSec < 72; iSec++){
 200         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 201         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 202         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 203         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 204
 205         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 206         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 207         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 208
 209         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 210         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 211         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 212         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 213
 214         if ( hQ != 0x0 ){
 215             TH2S *hNew = new TH2S(*hQ);
 216             hNew->SetDirectory(0);
 217             fHistoQArray.AddAt(hNew,iSec);
 218         }
 219         if ( hT0 != 0x0 ){
 220             TH2S *hNew = new TH2S(*hT0);
 221             hNew->SetDirectory(0);
 222             fHistoQArray.AddAt(hNew,iSec);
 223         }
 224         if ( hRMS != 0x0 ){
 225             TH2S *hNew = new TH2S(*hRMS);
 226             hNew->SetDirectory(0);
 227             fHistoQArray.AddAt(hNew,iSec);
 228         }
 229     }
 230     for (Int_t iEvent=0; iEvent<sig.fParamArrayEvent.GetSize(); iEvent++)
 231         fParamArrayEvent.AddAtAndExpand(sig.fParamArrayEvent.At(iEvent),iEvent);
 232     Int_t nrows = sig.fVEventTime.GetNrows();
 233     fVEventTime.ResizeTo(nrows);
 234     for (Int_t iEvent=0; iEvent<nrows; iEvent++)
 235         fVEventTime[iEvent] = sig.fVEventTime[iEvent];
 236
 237 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
 238 }
 239 //_____________________________________________________________________
 240 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
 241 {
 242   //
 243   // assignment operator
 244   //
 245   if (&source == this) return *this;
 246   new (this) AliTPCCalibCE(source);
 247
 248   return *this;
 249 }
 250 //_____________________________________________________________________
 251 AliTPCCalibCE::~AliTPCCalibCE()
 252 {
 253     //
 254     // destructor
 255     //
 256
 257     fCalRocArrayT0.Delete();
 258     fCalRocArrayQ.Delete();
 259     fCalRocArrayRMS.Delete();
 260
 261     fHistoQArray.Delete();
 262     fHistoT0Array.Delete();
 263     fHistoRMSArray.Delete();
 264
 265     fPadTimesArrayEvent.Delete();
 266     fPadQArrayEvent.Delete();
 267     fPadRMSArrayEvent.Delete();
 268     fPadPedestalArrayEvent.Delete();
 269
 270     if ( fDebugStreamer) delete fDebugStreamer;
 271 //    if ( fHTime0 ) delete fHTime0;
 272     delete fROC;
 273     delete fParam;
 274 }
 275 //_____________________________________________________________________
 276 Int_t AliTPCCalibCE::Update(const Int_t icsector, /*FOLD00*/
 277                                 const Int_t icRow,
 278                                 const Int_t icPad,
 279                                 const Int_t icTimeBin,
 280                                 const Float_t csignal)
 281 {
 282     //
 283     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 284     // no extra analysis necessary. Assumes knowledge of the signal shape!
 285     // assumes that it is looped over consecutive time bins of one pad
 286     //
 287     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 288
 289     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 290
 291     //init first pad and sector in this event
 292     if ( fCurrentChannel == -1 ) {
 293         fCurrentChannel = iChannel;
 294         fCurrentSector  = icsector;
 295         fCurrentRow     = icRow;
 296     }
 297
 298     //process last pad if we change to a new one
 299     if ( iChannel != fCurrentChannel ){
 300         ProcessPad();
 301         fCurrentChannel = iChannel;
 302         fCurrentSector  = icsector;
 303         fCurrentRow     = icRow;
 304     }
 305
 306     //fill signals for current pad
 307     fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
 308     if ( csignal > fMaxPadSignal ){
 309         fMaxPadSignal = csignal;
 310         fMaxTimeBin   = icTimeBin;
 311     }
 312     return 0;
 313 }
 314 //_____________________________________________________________________
 315 void AliTPCCalibCE::FindPedestal(Float_t part)
 316 {
 317     //
 318     // find pedestal and noise for the current pad. Use either database or
 319     // truncated mean with part*100%
 320     //
 321     Bool_t noPedestal = kTRUE;;
 322     if (fPedestalTPC&&fPadNoiseTPC){
 323         //use pedestal database
 324         //only load new pedestals if the sector has changed
 325         if ( fCurrentSector!=fLastSector ){
 326             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 327             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 328             fLastSector=fCurrentSector;
 329         }
 330
 331         if ( fPedestalROC&&fPadNoiseROC ){
 332             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
 333             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 334             noPedestal   = kFALSE;
 335         }
 336
 337     }
 338
 339     //if we are not running with pedestal database, or for the current sector there is no information
 340     //available, calculate the pedestal and noise on the fly
 341     if ( noPedestal ) {
 342         const Int_t kPedMax = 100;  //maximum pedestal value
 343         Float_t  max    =  0;
 344         Float_t  maxPos =  0;
 345         Int_t    median =  -1;
 346         Int_t    count0 =  0;
 347         Int_t    count1 =  0;
 348         //
 349         Float_t padSignal=0;
 350         //
 351         UShort_t histo[kPedMax];
 352         memset(histo,0,kPedMax*sizeof(UShort_t));
 353
 354         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; i++){
 355             padSignal = fPadSignal.GetMatrixArray()[i];
 356             if (padSignal<=0) continue;
 357             if (padSignal>max && i>10) {
 358                 max = padSignal;
 359                 maxPos = i;
 360             }
 361             if (padSignal>kPedMax-1) continue;
 362             histo[int(padSignal+0.5)]++;
 363             count0++;
 364         }
 365             //
 366         for (Int_t i=1; i<kPedMax; i++){
 367             if (count1<count0*0.5) median=i;
 368             count1+=histo[i];
 369         }
 370         // truncated mean
 371         //
 372         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 373         //
 374         for (Int_t idelta=1; idelta<10; idelta++){
 375             if (median-idelta<=0) continue;
 376             if (median+idelta>kPedMax) continue;
 377             if (count<part*count1){
 378                 count+=histo[median-idelta];
 379                 mean +=histo[median-idelta]*(median-idelta);
 380                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 381                 count+=histo[median+idelta];
 382                 mean +=histo[median+idelta]*(median+idelta);
 383                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 384             }
 385         }
 386         fPadPedestal = 0;
 387         fPadNoise    = 0;
 388         if ( count > 0 ) {
 389             mean/=count;
 390             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 391             fPadPedestal = mean;
 392             fPadNoise    = rms;
 393         }
 394     }
 395 }
 396 //_____________________________________________________________________
 397 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
 398 {
 399     //
 400     //  Find position, signal width and height of the CE signal (last signal)
 401     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 402     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 403     //
 404
 405     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 406     Int_t   cemaxpos       = 0;
 407     Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
 408     const Int_t    kCemin  = 4;             // range for the analysis of the ce signal +- channels from the peak
 409     const Int_t    kCemax  = 7;
 410
 411     Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
 412
 413     // find maximum closest to the sector mean from the last event
 414     for ( Int_t imax=0; imax<maxima.GetNrows(); imax++){
 415         Float_t tmean = (*((TVectorF*)(fTMeanArrayEvent[fTMeanArrayEvent.GetLast()])))[fCurrentSector];
 416             if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
 417                 minDist  = tmean-maxima[imax];
 418                 cemaxpos = (Int_t)maxima[imax];
 419             }
 420     }
 421
 422     if (cemaxpos!=0){
 423         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 424         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
 425             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 426             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
 427                 ceTime+=signal*(i+0.5);
 428                 ceRMS +=signal*(i+0.5)*(i+0.5);
 429                 ceQsum+=signal;
 430             }
 431         }
 432     }
 433     if (ceQmax&&ceQsum>ceSumThreshold) {
 434         ceTime/=ceQsum;
 435         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 436         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 437         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 438
 439         //Normalise Q to pad area of irocs
 440         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
 441
 442         ceQsum/=norm;
 443         fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
 444         fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
 445     } else {
 446         ceQmax=0;
 447         ceTime=0;
 448         ceRMS =0;
 449         ceQsum=0;
 450     }
 451     param[0] = ceQmax;
 452     param[1] = ceTime;
 453     param[2] = ceRMS;
 454     qSum     = ceQsum;
 455 }
 456 //_____________________________________________________________________
 457 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus)
 458 {
 459     //
 460     // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
 461     // and 'tplus' timebins after 'pos'
 462     //
 463     for (Int_t iTime = pos; iTime>pos-tminus; iTime--)
 464         if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
 465     for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; iTime++,iTime2++){
 466         if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
 467             iTime2++;
 468         if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
 469     }
 470     return kTRUE;
 471 }
 472 //_____________________________________________________________________
 473 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
 474 {
 475     //
 476     // Find local maxima on the pad signal and Histogram them
 477     //
 478     Float_t ceThreshold = 5.*fPadNoise;  // threshold for the signal
 479     Int_t   count       = 0;
 480     Int_t   tminus      = 2;
 481     Int_t   tplus       = 3;
 482     for (Int_t i=fLastTimeBin-tplus-1; i>=fFirstTimeBin+tminus; i--){
 483         if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,tminus,tplus) ){
 484             maxima.GetMatrixArray()[count++]=i;
 485             GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
 486         }
 487     }
 488 }
 489 //_____________________________________________________________________
 490 void AliTPCCalibCE::ProcessPad() /*FOLD00*/
 491 {
 492     //
 493     //  Process data of current pad
 494     //
 495     FindPedestal();
 496
 497     TVectorF maxima(10);
 498     FindLocalMaxima(maxima);
 499     if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
 500
 501
 502
 503     TVectorD param(3);
 504     Float_t  Qsum;
 505     FindCESignal(param, Qsum, maxima);
 506
 507     Double_t meanT  = param[1];
 508     Double_t sigmaT = param[2];
 509
 510     //Fill Event T0 counter
 511     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 512
 513     //Fill Q histogram
 514     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum), fCurrentChannel );
 515
 516     //Fill RMS histogram
 517     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 518
 519
 520     //Fill debugging info
 521     if ( fDebugLevel>0 ){
 522         (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
 523         (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
 524         (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=Qsum;
 525     }
 526
 527     ResetPad();
 528 }
 529 //_____________________________________________________________________
 530 void AliTPCCalibCE::EndEvent() /*FOLD00*/
 531 {
 532     //
 533     //  Process data of current pad
 534     //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
 535     //  before the EndEvent function to set the event timestamp and number!!!
 536     //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
 537     //  function was called
 538     //
 539
 540     //check if last pad has allready been processed, if not do so
 541     if ( fMaxTimeBin>-1 ) ProcessPad();
 542
 543     TVectorD param(3);
 544     TMatrixD dummy(3,3);
 545     TVectorF vMeanTime(72);
 546     TVectorF vMeanQ(72);
 547     AliTPCCalROC calIroc(0);
 548     AliTPCCalROC calOroc(36);
 549
 550     //find mean time0 offset for side A and C
 551     Double_t time0Side[2];       //time0 for side A:0 and C:0
 552     Double_t time0SideCount[2];  //time0 counter for side A:0 and C:0
 553     time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
 554     for ( Int_t iSec = 0; iSec<72; iSec++ ){
 555         time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
 556         time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
 557     }
 558     if ( time0SideCount[0] >0  )
 559         time0Side[0]/=time0SideCount[0];
 560     if ( time0SideCount[1] >0 )
 561         time0Side[1]/=time0SideCount[1];
 562     // end find time0 offset
 563
 564     //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
 565     for ( Int_t iSec = 0; iSec<72; iSec++ ){
 566
 567       //find median and then calculate the mean around it
 568         TH1S *hMeanT    = GetHistoTmean(iSec);
 569         if ( !hMeanT ) continue;
 570         Double_t entries = hMeanT->GetEntries();
 571         Double_t sum     = 0;
 572         Short_t *arr     = hMeanT->GetArray()+1;
 573         Int_t ibin=0;
 574         for ( ibin=0; ibin<hMeanT->GetNbinsX(); ibin++){
 575             sum+=arr[ibin];
 576             if ( sum>=(entries/2) ) break;
 577         }
 578         Int_t delta = 4;
 579         Int_t firstBin = fFirstTimeBin+ibin-delta;
 580         Int_t lastBin  = fFirstTimeBin+ibin+delta;
 581         if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
 582         if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
 583         Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
 584         vMeanTime.GetMatrixArray()[iSec]=median;
 585       // end find median
 586
 587         TVectorF *vTimes = GetPadTimesEvent(iSec);
 588         if ( !vTimes ) continue;
 589         AliTPCCalROC calIrocOutliers(0);
 590         AliTPCCalROC calOrocOutliers(36);
 591
 592         // calculate mean Q of the sector
 593         Float_t meanQ = 0;
 594         if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
 595         vMeanQ.GetMatrixArray()[iSec]=meanQ;
 596
 597         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++ ){
 598             Float_t Time  = (*vTimes).GetMatrixArray()[iChannel];
 599
 600             //set values for temporary roc calibration class
 601             if ( iSec < 36 ) {
 602                 calIroc.SetValue(iChannel, Time);
 603                 if ( Time == 0 ) calIrocOutliers.SetValue(iChannel,1);
 604
 605             } else {
 606                 calOroc.SetValue(iChannel, Time);
 607                 if ( Time == 0 ) calOrocOutliers.SetValue(iChannel,1);
 608             }
 609
 610             if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
 611                 GetHistoT0(iSec,kTRUE)->Fill( Time-time0Side[(iSec/18)%2],iChannel );
 612
 613
 614
 615             //-------------------------------  Debug start  ------------------------------
 616             if ( fDebugLevel>0 ){
 617                 if ( !fDebugStreamer ) {
 618                         //debug stream
 619                     TDirectory *backup = gDirectory;
 620                     fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
 621                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 622                 }
 623
 624                 Int_t row=0;
 625                 Int_t pad=0;
 626                 Int_t padc=0;
 627
 628                 Float_t Q   = (*GetPadQEvent(iSec))[iChannel];
 629                 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
 630
 631                 UInt_t channel=iChannel;
 632                 Int_t sector=iSec;
 633
 634                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 635                 pad = channel-fROC->GetRowIndexes(sector)[row];
 636                 padc = pad-(fROC->GetNPads(sector,row)/2);
 637
 638 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
 639 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
 640 //                                  fLastTimeBin-fFirstTimeBin,
 641 //                                  fFirstTimeBin,fLastTimeBin);
 642 //              h1->SetDirectory(0);
 643 //
 644 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
 645 //                  h1->Fill(i,fPadSignal(i));
 646
 647                 Double_t T0Sec = 0;
 648                 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
 649                     T0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
 650                 Double_t T0Side = time0Side[(iSec/18)%2];
 651                 (*fDebugStreamer) << "DataPad" <<
 652                     "Event=" << fNevents <<
 653                     "EventNumber=" << fRunNumber <<
 654                     "TimeStamp="   << fTimeStamp <<
 655                     "Sector="<< sector <<
 656                     "Row="   << row<<
 657                     "Pad="   << pad <<
 658                     "PadC="  << padc <<
 659                     "PadSec="<< channel <<
 660                     "Time0Sec="  << T0Sec <<
 661                     "Time0Side=" << T0Side <<
 662                     "Time="  << Time <<
 663                     "RMS="   << RMS <<
 664                     "Sum="   << Q <<
 665                     "MeanQ=" << meanQ <<
 666                     //              "hist.=" << h1 <<
 667                     "\n";
 668
 669                 //              delete h1;
 670
 671             }
 672             //-----------------------------  Debug end  ------------------------------
 673         }// end channel loop
 674         hMeanT->Reset();
 675
 676         TVectorD paramPol1(3);
 677         TVectorD paramPol2(6);
 678         TMatrixD matPol1(3,3);
 679         TMatrixD matPol2(6,6);
 680         Float_t  chi2Pol1=0;
 681         Float_t  chi2Pol2=0;
 682
 683         if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
 684             if ( iSec < 36 ){
 685                 calIroc.GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
 686                 calIroc.GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
 687             } else {
 688                 calOroc.GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
 689                 calOroc.GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
 690             }
 691
 692             GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
 693             GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
 694         }
 695 //      printf("events: %d -- size: %d\n",fNevents,GetParamArrayPol1(iSec)->GetSize());
 696
 697         //-------------------------------  Debug start  ------------------------------
 698         if ( fDebugLevel>0 ){
 699             if ( !fDebugStreamer ) {
 700                 //debug stream
 701                 TDirectory *backup = gDirectory;
 702                 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
 703                 if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 704             }
 705             (*fDebugStreamer) << "DataRoc" <<
 706                 "Event=" << fEvent <<
 707                 "EventNumber=" << fRunNumber <<
 708                 "TimeStamp="   << fTimeStamp <<
 709                 "Sector="<< iSec <<
 710                 "hMeanT.=" << hMeanT <<
 711                 "median=" << median <<
 712                 "paramPol1.=" << &paramPol1 <<
 713                 "paramPol2.=" << &paramPol2 <<
 714                 "matPol1.="   << &matPol1 <<
 715                 "matPol2.="   << &matPol2 <<
 716                 "chi2Pol1="   << chi2Pol1 <<
 717                 "chi2Pol2="   << chi2Pol2 <<
 718                 "\n";
 719         }
 720         //-------------------------------  Debug end  ------------------------------
 721     }// end sector loop
 722
 723 /*    AliMathBase::FitGaus(fHTime0->GetArray()+1,
 724                          fHTime0->GetNbinsX(),
 725                          fHTime0->GetXaxis()->GetXmin(),
 726                          fHTime0->GetXaxis()->GetXmax(),
 727                          &param, &dummy);*/
 728 //    fHTime0->Reset();
 729
 730     //    fParamArrayEvent.AddAtAndExpand(new TVectorD(param),fNevents);
 731     fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
 732     fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
 733     if ( fVEventTime.GetNrows() < fNevents ) {
 734         fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+1000));
 735         fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+1000));
 736     }
 737     fVEventTime[fNevents] = fTimeStamp;
 738     fVEventNumber[fNevents] = fRunNumber;
 739
 740     fNevents++;
 741     fOldRunNumber = fRunNumber;
 742
 743 }
 744 //_____________________________________________________________________
 745 Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
 746 {
 747   //
 748   // Event Processing loop - AliTPCRawStream
 749   // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
 750   //
 751
 752   rawStream->SetOldRCUFormat(fOldRCUformat);
 753
 754   ResetEvent();
 755
 756   Bool_t withInput = kFALSE;
 757
 758   while (rawStream->Next()) {
 759
 760       Int_t isector  = rawStream->GetSector();                       //  current sector
 761       Int_t iRow     = rawStream->GetRow();                          //  current row
 762       Int_t iPad     = rawStream->GetPad();                          //  current pad
 763       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 764       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 765
 766       Update(isector,iRow,iPad,iTimeBin,signal);
 767       withInput = kTRUE;
 768   }
 769
 770   if (withInput){
 771       EndEvent();
 772   }
 773
 774   return withInput;
 775 }
 776 //_____________________________________________________________________
 777 Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
 778 {
 779   //
 780   //  Event processing loop - AliRawReader
 781   //
 782
 783
 784     AliTPCRawStream rawStream(rawReader);
 785     AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
 786     if (eventHeader){
 787         fTimeStamp   = eventHeader->Get("Timestamp");
 788         fRunNumber = eventHeader->Get("RunNb");
 789     }
 790
 791
 792   rawReader->Select("TPC");
 793
 794   return ProcessEvent(&rawStream);
 795 }
 796 //_____________________________________________________________________
 797 Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
 798 {
 799   //
 800   //  Event processing loop - date event
 801   //
 802     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 803     Bool_t result=ProcessEvent(rawReader);
 804     delete rawReader;
 805     return result;
 806
 807 }
 808 //_____________________________________________________________________
 809 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
 810                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 811                                   Char_t *type, Bool_t force)
 812 {
 813     //
 814     // return pointer to TH2S histogram of 'type'
 815     // if force is true create a new histogram if it doesn't exist allready
 816     //
 817     if ( !force || arr->UncheckedAt(sector) )
 818         return (TH2S*)arr->UncheckedAt(sector);
 819
 820     // if we are forced and histogram doesn't yes exist create it
 821     Char_t name[255], title[255];
 822
 823     sprintf(name,"hCalib%s%.2d",type,sector);
 824     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 825
 826     // new histogram with Q calib information. One value for each pad!
 827     TH2S* hist = new TH2S(name,title,
 828                           nbinsY, ymin, ymax,
 829                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 830     hist->SetDirectory(0);
 831     arr->AddAt(hist,sector);
 832     return hist;
 833 }
 834 //_____________________________________________________________________
 835 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
 836 {
 837     //
 838     // return pointer to T0 histogram
 839     // if force is true create a new histogram if it doesn't exist allready
 840     //
 841     TObjArray *arr = &fHistoT0Array;
 842     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 843 }
 844 //_____________________________________________________________________
 845 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
 846 {
 847     //
 848     // return pointer to Q histogram
 849     // if force is true create a new histogram if it doesn't exist allready
 850     //
 851     TObjArray *arr = &fHistoQArray;
 852     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 853 }
 854 //_____________________________________________________________________
 855 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
 856 {
 857     //
 858     // return pointer to Q histogram
 859     // if force is true create a new histogram if it doesn't exist allready
 860     //
 861     TObjArray *arr = &fHistoRMSArray;
 862     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 863 }
 864 //_____________________________________________________________________
 865 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
 866                               Char_t *type, Bool_t force)
 867 {
 868     //
 869     // return pointer to TH1S histogram
 870     // if force is true create a new histogram if it doesn't exist allready
 871     //
 872     if ( !force || arr->UncheckedAt(sector) )
 873         return (TH1S*)arr->UncheckedAt(sector);
 874
 875     // if we are forced and histogram doesn't yes exist create it
 876     Char_t name[255], title[255];
 877
 878     sprintf(name,"hCalib%s%.2d",type,sector);
 879     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 880
 881     // new histogram with Q calib information. One value for each pad!
 882     TH1S* hist = new TH1S(name,title,
 883                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
 884     hist->SetDirectory(0);
 885     arr->AddAt(hist,sector);
 886     return hist;
 887 }
 888 //_____________________________________________________________________
 889 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
 890 {
 891     //
 892     // return pointer to Q histogram
 893     // if force is true create a new histogram if it doesn't exist allready
 894     //
 895     TObjArray *arr = &fHistoTmean;
 896     return GetHisto(sector, arr, "LastTmean", force);
 897 }
 898 //_____________________________________________________________________
 899 TVectorF* AliTPCCalibCE::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
 900 {
 901     //
 902     // return pointer to Pad Info from 'arr' for the current event and sector
 903     // if force is true create it if it doesn't exist allready
 904     //
 905     if ( !force || arr->UncheckedAt(sector) )
 906         return (TVectorF*)arr->UncheckedAt(sector);
 907
 908     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 909     arr->AddAt(vect,sector);
 910     return vect;
 911 }
 912 //_____________________________________________________________________
 913 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
 914 {
 915     //
 916     // return pointer to Pad Times Array for the current event and sector
 917     // if force is true create it if it doesn't exist allready
 918     //
 919     TObjArray *arr = &fPadTimesArrayEvent;
 920     return GetPadInfoEvent(sector,arr,force);
 921 }
 922 //_____________________________________________________________________
 923 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
 924 {
 925     //
 926     // return pointer to Pad Q Array for the current event and sector
 927     // if force is true create it if it doesn't exist allready
 928     // for debugging purposes only
 929     //
 930
 931     TObjArray *arr = &fPadQArrayEvent;
 932     return GetPadInfoEvent(sector,arr,force);
 933 }
 934 //_____________________________________________________________________
 935 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
 936 {
 937     //
 938     // return pointer to Pad RMS Array for the current event and sector
 939     // if force is true create it if it doesn't exist allready
 940     // for debugging purposes only
 941     //
 942     TObjArray *arr = &fPadRMSArrayEvent;
 943     return GetPadInfoEvent(sector,arr,force);
 944 }
 945 //_____________________________________________________________________
 946 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
 947 {
 948     //
 949     // return pointer to Pad RMS Array for the current event and sector
 950     // if force is true create it if it doesn't exist allready
 951     // for debugging purposes only
 952     //
 953     TObjArray *arr = &fPadPedestalArrayEvent;
 954     return GetPadInfoEvent(sector,arr,force);
 955 }
 956 //_____________________________________________________________________
 957 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
 958 {
 959     //
 960     // return pointer to ROC Calibration
 961     // if force is true create a new histogram if it doesn't exist allready
 962     //
 963     if ( !force || arr->UncheckedAt(sector) )
 964         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 965
 966     // if we are forced and histogram doesn't yes exist create it
 967
 968     // new AliTPCCalROC for T0 information. One value for each pad!
 969     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 970     //init values
 971     for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
 972         croc->SetValue(iChannel, 0);
 973     }
 974     arr->AddAt(croc,sector);
 975     return croc;
 976 }
 977 //_____________________________________________________________________
 978 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
 979 {
 980     //
 981     // return pointer to Carge ROC Calibration
 982     // if force is true create a new histogram if it doesn't exist allready
 983     //
 984     TObjArray *arr = &fCalRocArrayT0;
 985     return GetCalRoc(sector, arr, force);
 986 }
 987 //_____________________________________________________________________
 988 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
 989 {
 990     //
 991     // return pointer to T0 ROC Calibration
 992     // if force is true create a new histogram if it doesn't exist allready
 993     //
 994     TObjArray *arr = &fCalRocArrayQ;
 995     return GetCalRoc(sector, arr, force);
 996 }
 997 //_____________________________________________________________________
 998 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
 999 {
1000     //
1001     // return pointer to signal width ROC Calibration
1002     // if force is true create a new histogram if it doesn't exist allready
1003     //
1004     TObjArray *arr = &fCalRocArrayRMS;
1005     return GetCalRoc(sector, arr, force);
1006 }
1007 //_____________________________________________________________________
1008 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1009 {
1010     //
1011     // return pointer to Outliers
1012     // if force is true create a new histogram if it doesn't exist allready
1013     //
1014     TObjArray *arr = &fCalRocArrayOutliers;
1015     return GetCalRoc(sector, arr, force);
1016 }
1017 //_____________________________________________________________________
1018 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force)
1019 {
1020     //
1021     // return pointer to TObjArray of fit parameters
1022     // if force is true create a new histogram if it doesn't exist allready
1023     //
1024     if ( !force || arr->UncheckedAt(sector) )
1025         return (TObjArray*)arr->UncheckedAt(sector);
1026
1027     // if we are forced and array doesn't yes exist create it
1028
1029     // new TObjArray for parameters
1030     TObjArray *newArr = new TObjArray;
1031     arr->AddAt(newArr,sector);
1032     return newArr;
1033 }
1034 //_____________________________________________________________________
1035 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1036 {
1037     //
1038     // return pointer to TObjArray of fit parameters from plane fit
1039     // if force is true create a new histogram if it doesn't exist allready
1040     //
1041     TObjArray *arr = &fParamArrayEventPol1;
1042     return GetParamArray(sector, arr, force);
1043 }
1044 //_____________________________________________________________________
1045 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1046 {
1047     //
1048     // return pointer to TObjArray of fit parameters from parabola fit
1049     // if force is true create a new histogram if it doesn't exist allready
1050     //
1051     TObjArray *arr = &fParamArrayEventPol2;
1052     return GetParamArray(sector, arr, force);
1053 }
1054 //_____________________________________________________________________
1055 void AliTPCCalibCE::ResetEvent() /*FOLD00*/
1056 {
1057     //
1058     //  Reset global counters  -- Should be called before each event is processed
1059     //
1060     fLastSector=-1;
1061     fCurrentSector=-1;
1062     fCurrentRow=-1;
1063     fCurrentChannel=-1;
1064
1065     ResetPad();
1066
1067     fPadTimesArrayEvent.Delete();
1068     fPadQArrayEvent.Delete();
1069     fPadRMSArrayEvent.Delete();
1070     fPadPedestalArrayEvent.Delete();
1071
1072     for ( Int_t i=0; i<72; i++ ){
1073         fVTime0Offset.GetMatrixArray()[i]=0;
1074         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1075         fVMeanQ.GetMatrixArray()[i]=0;
1076         fVMeanQCounter.GetMatrixArray()[i]=0;
1077     }
1078 }
1079 //_____________________________________________________________________
1080 void AliTPCCalibCE::ResetPad() /*FOLD00*/
1081 {
1082     //
1083     //  Reset pad infos -- Should be called after a pad has been processed
1084     //
1085     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
1086         fPadSignal.GetMatrixArray()[i] = 0;
1087     fMaxTimeBin   = -1;
1088     fMaxPadSignal = -1;
1089     fPadPedestal  = -1;
1090     fPadNoise     = -1;
1091 }
1092 //_____________________________________________________________________
1093 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter) /*FOLD00*/
1094 {
1095     //
1096     // Make graph from fit parameters of pol1 or pol2 fit
1097     // xVariable:    0-run time, 1-run number, 2-internal event counter
1098     // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 2-mean Q
1099     // fitParameter: fit parameter ( 0-2 for pol1, 0-5 for pol2, 0 for mean time )
1100     //
1101
1102     Double_t *x = new Double_t[fNevents];
1103     Double_t *y = new Double_t[fNevents];
1104
1105     TVectorD *xVar = 0x0;
1106     TObjArray *aType = 0x0;
1107     Int_t npoints=0;
1108
1109     // sanity checks
1110     if ( (sector<0) || (sector>71) )      return 0x0;
1111     if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1112     if ( (fitType<0) || (fitType>3) )     return 0x0;
1113     if ( fitType==0 ){
1114         if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1115         aType = &fParamArrayEventPol1;
1116         if ( aType->At(sector)==0x0 ) return 0x0;
1117     }
1118     else if ( fitType==1 ){
1119         if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1120         aType = &fParamArrayEventPol2;
1121         if ( aType->At(sector)==0x0 ) return 0x0;
1122     }
1123
1124
1125     if ( xVariable == 0 ) xVar = &fVEventTime;
1126     if ( xVariable == 1 ) xVar = &fVEventNumber;
1127     if ( xVariable == 2 ) {
1128         xVar = new TVectorD(fNevents);
1129         for ( Int_t i=0;i<fNevents; i++) (*xVar)[i]=i;
1130     }
1131
1132     for (Int_t ievent =0; ievent<fNevents; ievent++){
1133         if ( fitType<2 ){
1134             TObjArray *events = (TObjArray*)(aType->At(sector));
1135             if ( events->GetSize()<=ievent ) break;
1136             TVectorD *v = (TVectorD*)(events->At(ievent));
1137             if ( v!=0x0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1138         } else if (fitType == 2) {
1139             Double_t yValue=(*((TVectorF*)(fTMeanArrayEvent[ievent])))[sector];
1140             if ( yValue>0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=yValue;npoints++;}
1141         }else if (fitType == 3) {
1142             Double_t yValue=(*((TVectorF*)(fQMeanArrayEvent[ievent])))[sector];
1143             if ( yValue>0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=yValue;npoints++;}
1144         }
1145     }
1146
1147     TGraph *gr = new TGraph(npoints);
1148     //sort xVariable increasing
1149     Int_t    *sortIndex = new Int_t[npoints];
1150     TMath::Sort(npoints,x,sortIndex);
1151     for (Int_t i=0;i<npoints;i++){
1152         gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1153     }
1154
1155
1156     if ( xVariable == 2 ) delete xVar;
1157     delete x;
1158     delete y;
1159     delete sortIndex;
1160     return gr;
1161 }
1162 //_____________________________________________________________________
1163 void AliTPCCalibCE::Analyse()
1164 {
1165     //
1166     //  Calculate calibration constants
1167     //
1168
1169     TVectorD paramQ(3);
1170     TVectorD paramT0(3);
1171     TVectorD paramRMS(3);
1172     TMatrixD dummy(3,3);
1173
1174     for (Int_t iSec=0; iSec<72; iSec++){
1175         TH2S *hT0 = GetHistoT0(iSec);
1176         if (!hT0 ) continue;
1177
1178         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1179         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1180         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1181         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1182
1183         TH2S *hQ   = GetHistoQ(iSec);
1184         TH2S *hRMS = GetHistoRMS(iSec);
1185
1186         Short_t *array_hQ   = hQ->GetArray();
1187         Short_t *array_hT0  = hT0->GetArray();
1188         Short_t *array_hRMS = hRMS->GetArray();
1189
1190         UInt_t nChannels = fROC->GetNChannels(iSec);
1191
1192         //debug
1193         Int_t row=0;
1194         Int_t pad=0;
1195         Int_t padc=0;
1196         //! debug
1197
1198         for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
1199
1200
1201             Float_t cogTime0 = -1000;
1202             Float_t cogQ     = -1000;
1203             Float_t cogRMS   = -1000;
1204             Float_t cogOut   = 0;
1205
1206
1207             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1208             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1209             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1210
1211 /*
1212             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1213             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1214             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1215             cogQ     = paramQ[1];
1216             cogTime0 = paramT0[1];
1217             cogRMS   = paramRMS[1];
1218 */
1219             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1220             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1221             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1222
1223
1224
1225             /*
1226             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1227                 cogOut = 1;
1228                 cogTime0 = 0;
1229                 cogQ     = 0;
1230                 cogRMS   = 0;
1231             }
1232 */
1233             rocQ->SetValue(iChannel, cogQ*cogQ);
1234             rocT0->SetValue(iChannel, cogTime0);
1235             rocRMS->SetValue(iChannel, cogRMS);
1236             rocOut->SetValue(iChannel, cogOut);
1237
1238
1239             //debug
1240             if ( fDebugLevel > 0 ){
1241                 if ( !fDebugStreamer ) {
1242                         //debug stream
1243                     TDirectory *backup = gDirectory;
1244                     fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1245                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1246                 }
1247
1248                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1249                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1250                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1251
1252                 (*fDebugStreamer) << "DataEnd" <<
1253                     "Sector="  << iSec      <<
1254                     "Pad="     << pad       <<
1255                     "PadC="    << padc      <<
1256                     "Row="     << row       <<
1257                     "PadSec="  << iChannel   <<
1258                     "Q="       << cogQ      <<
1259                     "T0="      << cogTime0  <<
1260                     "RMS="     << cogRMS    <<
1261                     "\n";
1262             }
1263             //! debug
1264
1265         }
1266
1267     }
1268     fDebugStreamer->GetFile()->Write();
1269 //    delete fDebugStreamer;
1270 //    fDebugStreamer = 0x0;
1271 }
1272 //_____________________________________________________________________
1273 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1274 {
1275     //
1276     //  Write class to file
1277     //
1278
1279     TString sDir(dir);
1280     TString option;
1281
1282     if ( append )
1283         option = "update";
1284     else
1285         option = "recreate";
1286
1287     TDirectory *backup = gDirectory;
1288     TFile f(filename,option.Data());
1289     f.cd();
1290     if ( !sDir.IsNull() ){
1291         f.mkdir(sDir.Data());
1292         f.cd(sDir);
1293     }
1294     this->Write();
1295     f.Close();
1296
1297     if ( backup ) backup->cd();
1298 }