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 /* $Id$ */
  17
  18 ////////////////////////////////////////////////////////////////////////////////////////
  19 //                                                                                    //
  20 //             Implementation of the TPC Central Electrode calibration                //
  21 //                                                                                    //
  22 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch  //
  23 //                                                                                    //
  24 ////////////////////////////////////////////////////////////////////////////////////////
  25 //
  26 //
  27 // *************************************************************************************
  28 // *                                Class Description                                  *
  29 // *************************************************************************************
  30 //
  31 /* BEGIN_HTML
  32  <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
  33  using laser runs.</h4>
  34
  35  The information retrieved is
  36  <ul style="list-style-type: square;">
  37    <li>Time arrival from the CE</li>
  38    <li>Signal width</li>
  39    <li>Signal sum</li>
  40  </ul>
  41
  42 <h4>Overview:</h4>
  43  <ol style="list-style-type: upper-roman;">
  44    <li><a href="#working">Working principle</a></li>
  45    <li><a href="#user">User interface for filling data</a></li>
  46    <li><a href="#info">Stored information</a></li>
  47  </ol>
  48
  49  <h3><a name="working">I. Working principle</a></h3>
  50
  51  <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
  52  (see below). These in the end call the Update(...) function.</h4>
  53
  54  <ul style="list-style-type: square;">
  55    <li>the Update(...) function:<br />
  56        In this function the array fPadSignal is filled with the adc signals between the specified range
  57        fFirstTimeBin and fLastTimeBin for the current pad.
  58        before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
  59        stored in fPadSignal.
  60    </li>
  61    <ul style="list-style-type: square;">
  62    <li>the ProcessPad() function:</li>
  63    <ol style="list-style-type: decimal;">
  64      <li>Find Pedestal and Noise information</li>
  65      <ul style="list-style-type: square;">
  66        <li>use database information which has to be set by calling<br />
  67            SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
  68        <li>if no information from the pedestal data base
  69            is available the informaion is calculated on the fly
  70            ( see FindPedestal() function )</li>
  71      </ul>
  72      <li>Find local maxima of the pad signal</li>
  73      <ul style="list-style-type: square;">
  74        <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
  75        have been observed ( see FindLocalMaxima(...) )</li>
  76      </ul>
  77      <li>Find the CE signal information</li>
  78      <ul style="list-style-type: square;">
  79        <li>to find the position of the CE signal the Tmean information from the previos event is used
  80            as the CE signal the local maximum closest to this Tmean is identified</li>
  81        <li>calculate  mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
  82            the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
  83      </ul>
  84      <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
  85      <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
  86      </ol>
  87    </ul>
  88  </ul>
  89
  90  <h4>At the end of each event the EndEvent() function is called</h4>
  91
  92  <ul style="list-style-type: square;">
  93    <li>the EndEvent() function:</li>
  94    <ul style="list-style-type: square;">
  95      <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
  96          This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
  97      <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
  98      <li>calculate Mean Q</li>
  99      <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
 100    </ul>
 101  </ul>
 102
 103  <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
 104   <ul style="list-style-type: square;">
 105   <li>the Analyse() function:</li>
 106     <ul style="list-style-type: square;">
 107       <li>calculate the mean values of T0, RMS, Q for each pad, using
 108           the AliMathBase::GetCOG(...) function</li>
 109       <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
 110          (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
 111     </ul>
 112   </ul>
 113
 114  <h3><a name="user">II. User interface for filling data</a></h3>
 115
 116  <h4>To Fill information one of the following functions can be used:</h4>
 117
 118  <ul style="list-style-type: none;">
 119   <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
 120     <ul style="list-style-type: square;">
 121       <li>process Date event</li>
 122       <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
 123     </ul>
 124     <br />
 125
 126   <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
 127     <ul style="list-style-type: square;">
 128       <li>process AliRawReader event</li>
 129       <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
 130     </ul>
 131     <br />
 132
 133   <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
 134     <ul style="list-style-type: square;">
 135       <li>process event from AliTPCRawStream</li>
 136       <li>call Update function for signal filling</li>
 137     </ul>
 138     <br />
 139
 140   <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
 141               iPad,  const Int_t iTimeBin, const Float_t signal);</li>
 142     <ul style="list-style-type: square;">
 143       <li>directly  fill signal information (sector, row, pad, time bin, pad)
 144           to the reference histograms</li>
 145     </ul>
 146  </ul>
 147
 148  <h4>It is also possible to merge two independently taken calibrations using the function</h4>
 149
 150  <ul style="list-style-type: none;">
 151   <li> void Merge(AliTPCCalibSignal *sig)</li>
 152     <ul style="list-style-type: square;">
 153       <li>copy histograms in 'sig' if they do not exist in this instance</li>
 154       <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
 155       <li>After merging call Analyse again!</li>
 156     </ul>
 157  </ul>
 158
 159
 160  <h4>example: filling data using root raw data:</h4>
 161  <pre>
 162  void fillCE(Char_t *filename)
 163  {
 164     rawReader = new AliRawReaderRoot(fileName);
 165     if ( !rawReader ) return;
 166     AliTPCCalibCE *calib = new AliTPCCalibCE;
 167     while (rawReader->NextEvent()){
 168       calib->ProcessEvent(rawReader);
 169     }
 170     calib->Analyse();
 171     calib->DumpToFile("CEData.root");
 172     delete rawReader;
 173     delete calib;
 174  }
 175  </pre>
 176
 177  <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
 178
 179  <h4><a name="info:stored">III.1 Stored information</a></h4>
 180  <ul style="list-style-type: none;">
 181   <li>Histograms:</li>
 182   <ul style="list-style-type: none;">
 183     <li>For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum])
 184         is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
 185         stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
 186   </ul>
 187   <br />
 188
 189  <li>Calibration Data:</li>
 190  <ul style="list-style-type: none;">
 191       <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
 192           the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
 193           fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
 194           is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
 195  </ul>
 196  <br />
 197
 198  <li>For each event the following information is stored:</li>
 199
 200  <ul style="list-style-type: square;">
 201    <li>event time ( TVectorD  fVEventTime )</li>
 202    <li>event id   ( TVectorD  fVEventNumber )</li>
 203    <br />
 204    <li>mean arrival time for each ROC                ( TObjArray fTMeanArrayEvent )</li>
 205    <li>mean Q for each ROC                           ( TObjArray fQMeanArrayEvent )</li>
 206    <li>parameters of a plane fit for each ROC        ( TObjArray fParamArrayEventPol1 )</li>
 207    <li>parameters of a 2D parabola fit for each ROC  ( TObjArray fParamArrayEventPol2 )</li>
 208   </ul>
 209  </ul>
 210
 211  <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
 212  <ul style="list-style-type: none;">
 213   <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
 214     <ul style="list-style-type: square;">
 215       <li>TH2F *GetHistoT0(Int_t sector);</li>
 216       <li>TH2F *GetHistoRMS(Int_t sector);</li>
 217       <li>TH2F *GetHistoQ(Int_t sector);</li>
 218     </ul>
 219     <br />
 220
 221   <li>Accessing the calibration storage objects:</li>
 222     <ul style="list-style-type: square;">
 223       <li>AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values</li>
 224       <li>AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values</li>
 225       <li>AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values</li>
 226     </ul>
 227     <br />
 228
 229   <li>Accessin the event by event information:</li>
 230     <ul style="list-style-type: square;">
 231       <li>The event by event information can be displayed using the</li>
 232       <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
 233       <li>which creates a graph from the specified variables</li>
 234     </ul>
 235   </ul>
 236
 237   <h4>example for visualisation:</h4>
 238   <pre>
 239   //if the file "CEData.root" was created using the above example one could do the following:
 240   TFile fileCE("CEData.root")
 241   AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
 242   ce->GetCalRocT0(0)->Draw("colz");
 243   ce->GetCalRocRMS(0)->Draw("colz");
 244
 245   //or use the AliTPCCalPad functionality:
 246   AliTPCCalPad padT0(ped->GetCalPadT0());
 247   AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
 248   padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
 249   padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
 250
 251   //display event by event information:
 252   //Draw mean arrival time as a function of the event time for oroc sector A00
 253   ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
 254   //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
 255   ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");
 256   </pre>
 257 END_HTML */
 258 //////////////////////////////////////////////////////////////////////////////////////
 259
 260
 261 //Root includes
 262 #include <TObjArray.h>
 263 #include <TH1F.h>
 264 #include <TH2S.h>
 265 #include <TString.h>
 266 #include <TVectorF.h>
 267 #include <TVectorD.h>
 268 #include <TMatrixD.h>
 269 #include <TMath.h>
 270 #include <TGraph.h>
 271 #include <TString.h>
 272 #include <TMap.h>
 273 #include <TDirectory.h>
 274 #include <TSystem.h>
 275 #include <TFile.h>
 276
 277 //AliRoot includes
 278 #include "AliLog.h"
 279 #include "AliRawReader.h"
 280 #include "AliRawReaderRoot.h"
 281 #include "AliRawReaderDate.h"
 282 #include "AliRawEventHeaderBase.h"
 283 #include "AliTPCRawStream.h"
 284 #include "AliTPCRawStreamFast.h"
 285 #include "AliTPCcalibDB.h"
 286 #include "AliTPCCalROC.h"
 287 #include "AliTPCCalPad.h"
 288 #include "AliTPCROC.h"
 289 #include "AliTPCParam.h"
 290 #include "AliTPCCalibCE.h"
 291 #include "AliMathBase.h"
 292 #include "TTreeStream.h"
 293
 294 //date
 295 #include "event.h"
 296 ClassImp(AliTPCCalibCE)
 297
 298
 299 AliTPCCalibCE::AliTPCCalibCE() :
 300     TObject(),
 301     fFirstTimeBin(650),
 302     fLastTimeBin(1000),
 303     fNbinsT0(200),
 304     fXminT0(-5),
 305     fXmaxT0(5),
 306     fNbinsQ(200),
 307     fXminQ(1),
 308     fXmaxQ(40),
 309     fNbinsRMS(100),
 310     fXminRMS(0.1),
 311     fXmaxRMS(5.1),
 312     fPeakMinus(2),
 313     fPeakPlus(3),
 314     fNoiseThresholdMax(5.),
 315     fNoiseThresholdSum(8.),
 316     fIsZeroSuppressed(kFALSE),
 317     fLastSector(-1),
 318     fSecRejectRatio(.4),
 319     fROC(AliTPCROC::Instance()),
 320     fMapping(NULL),
 321     fParam(new AliTPCParam),
 322     fPedestalTPC(0x0),
 323     fPadNoiseTPC(0x0),
 324     fPedestalROC(0x0),
 325     fPadNoiseROC(0x0),
 326     fCalRocArrayT0(72),
 327     fCalRocArrayT0Err(72),
 328     fCalRocArrayQ(72),
 329     fCalRocArrayRMS(72),
 330     fCalRocArrayOutliers(72),
 331     fHistoQArray(72),
 332     fHistoT0Array(72),
 333     fHistoRMSArray(72),
 334     fMeanT0rms(0),
 335     fMeanQrms(0),
 336     fMeanRMSrms(0),
 337     fHistoTmean(72),
 338     fParamArrayEventPol1(72),
 339     fParamArrayEventPol2(72),
 340     fTMeanArrayEvent(72),
 341     fQMeanArrayEvent(72),
 342     fVEventTime(10),
 343     fVEventNumber(10),
 344     fNevents(0),
 345     fTimeStamp(0),
 346     fEventId(-1),
 347     fRunNumber(-1),
 348     fOldRunNumber(-1),
 349     fPadTimesArrayEvent(72),
 350     fPadQArrayEvent(72),
 351     fPadRMSArrayEvent(72),
 352     fPadPedestalArrayEvent(72),
 353     fCurrentChannel(-1),
 354     fCurrentSector(-1),
 355     fCurrentRow(-1),
 356     fMaxPadSignal(-1),
 357     fMaxTimeBin(-1),
 358     fPadSignal(1024),
 359     fPadPedestal(0),
 360     fPadNoise(0),
 361     fVTime0Offset(72),
 362     fVTime0OffsetCounter(72),
 363     fVMeanQ(72),
 364     fVMeanQCounter(72),
 365     fCurrentCETimeRef(0),
 366 //    fEvent(-1),
 367     fDebugStreamer(0x0),
 368     fDebugLevel(0)
 369 {
 370     //
 371     // AliTPCSignal default constructor
 372     //
 373 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
 374     fParam->Update();
 375 }
 376 //_____________________________________________________________________
 377 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
 378     TObject(sig),
 379     fFirstTimeBin(sig.fFirstTimeBin),
 380     fLastTimeBin(sig.fLastTimeBin),
 381     fNbinsT0(sig.fNbinsT0),
 382     fXminT0(sig.fXminT0),
 383     fXmaxT0(sig.fXmaxT0),
 384     fNbinsQ(sig.fNbinsQ),
 385     fXminQ(sig.fXminQ),
 386     fXmaxQ(sig.fXmaxQ),
 387     fNbinsRMS(sig.fNbinsRMS),
 388     fXminRMS(sig.fXminRMS),
 389     fXmaxRMS(sig.fXmaxRMS),
 390     fPeakMinus(sig.fPeakMinus),
 391     fPeakPlus(sig.fPeakPlus),
 392     fNoiseThresholdMax(sig.fNoiseThresholdMax),
 393     fNoiseThresholdSum(sig.fNoiseThresholdSum),
 394     fIsZeroSuppressed(sig.fIsZeroSuppressed),
 395     fLastSector(-1),
 396     fSecRejectRatio(.4),
 397     fROC(AliTPCROC::Instance()),
 398     fMapping(NULL),
 399     fParam(new AliTPCParam),
 400     fPedestalTPC(0x0),
 401     fPadNoiseTPC(0x0),
 402     fPedestalROC(0x0),
 403     fPadNoiseROC(0x0),
 404     fCalRocArrayT0(72),
 405     fCalRocArrayT0Err(72),
 406     fCalRocArrayQ(72),
 407     fCalRocArrayRMS(72),
 408     fCalRocArrayOutliers(72),
 409     fHistoQArray(72),
 410     fHistoT0Array(72),
 411     fHistoRMSArray(72),
 412     fMeanT0rms(sig.fMeanT0rms),
 413     fMeanQrms(sig.fMeanQrms),
 414     fMeanRMSrms(sig.fMeanRMSrms),
 415     fHistoTmean(72),
 416     fParamArrayEventPol1(72),
 417     fParamArrayEventPol2(72),
 418     fTMeanArrayEvent(72),
 419     fQMeanArrayEvent(72),
 420     fVEventTime(1000),
 421     fVEventNumber(1000),
 422     fNevents(sig.fNevents),
 423     fTimeStamp(0),
 424     fEventId(-1),
 425     fRunNumber(-1),
 426     fOldRunNumber(-1),
 427     fPadTimesArrayEvent(72),
 428     fPadQArrayEvent(72),
 429     fPadRMSArrayEvent(72),
 430     fPadPedestalArrayEvent(72),
 431     fCurrentChannel(-1),
 432     fCurrentSector(-1),
 433     fCurrentRow(-1),
 434     fMaxPadSignal(-1),
 435     fMaxTimeBin(-1),
 436     fPadSignal(1024),
 437     fPadPedestal(0),
 438     fPadNoise(0),
 439     fVTime0Offset(72),
 440     fVTime0OffsetCounter(72),
 441     fVMeanQ(72),
 442     fVMeanQCounter(72),
 443     fCurrentCETimeRef(0),
 444 //    fEvent(-1),
 445     fDebugStreamer(0x0),
 446     fDebugLevel(sig.fDebugLevel)
 447 {
 448   //
 449     // AliTPCSignal copy constructor
 450   //
 451
 452   for (Int_t iSec = 0; iSec < 72; ++iSec){
 453     const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 454     const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 455     const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 456     const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 457
 458     const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 459     const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 460     const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 461
 462     if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 463     if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 464     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 465     if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 466
 467     if ( hQ != 0x0 ){
 468 //          TDirectory *dir = hQ->GetDirectory();
 469 //          hQ->SetDirectory(0);
 470       TH2S *hNew = new TH2S(*hQ);
 471       hNew->SetDirectory(0);
 472       fHistoQArray.AddAt(hNew,iSec);
 473 //            hQ->SetDirectory(dir);
 474     }
 475     if ( hT0 != 0x0 ){
 476 //          TDirectory *dir = hT0->GetDirectory();
 477 //          hT0->SetDirectory(0);
 478       TH2S *hNew = new TH2S(*hT0);
 479       hNew->SetDirectory(0);
 480       fHistoT0Array.AddAt(hNew,iSec);
 481 //            hT0->SetDirectory(dir);
 482     }
 483     if ( hRMS != 0x0 ){
 484 //          TDirectory *dir = hRMS->GetDirectory();
 485 //          hRMS->SetDirectory(0);
 486       TH2S *hNew = new TH2S(*hRMS);
 487       hNew->SetDirectory(0);
 488       fHistoRMSArray.AddAt(hNew,iSec);
 489 //            hRMS->SetDirectory(dir);
 490     }
 491   }
 492
 493     //copy fit parameters event by event
 494   TObjArray *arr=0x0;
 495   for (Int_t iSec=0; iSec<72; ++iSec){
 496     arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
 497     if ( arr ){
 498       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 499       fParamArrayEventPol1.AddAt(arrEvents, iSec);
 500       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 501         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 502           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 503     }
 504
 505     arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
 506     if ( arr ){
 507       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 508       fParamArrayEventPol2.AddAt(arrEvents, iSec);
 509       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 510         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 511           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 512     }
 513
 514     TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
 515     TVectorF *vMeanQ    = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
 516     if ( vMeanTime )
 517       fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
 518     if ( vMeanQ )
 519       fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
 520   }
 521
 522
 523   fVEventTime.ResizeTo(sig.fVEventTime);
 524   fVEventNumber.ResizeTo(sig.fVEventNumber);
 525   fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
 526   fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
 527
 528   fParam->Update();
 529 }
 530 //_____________________________________________________________________
 531 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
 532   TObject(),
 533   fFirstTimeBin(650),
 534   fLastTimeBin(1000),
 535   fNbinsT0(200),
 536   fXminT0(-5),
 537   fXmaxT0(5),
 538   fNbinsQ(200),
 539   fXminQ(1),
 540   fXmaxQ(40),
 541   fNbinsRMS(100),
 542   fXminRMS(0.1),
 543   fXmaxRMS(5.1),
 544   fPeakMinus(2),
 545   fPeakPlus(3),
 546   fNoiseThresholdMax(5.),
 547   fNoiseThresholdSum(8.),
 548   fIsZeroSuppressed(kFALSE),
 549   fLastSector(-1),
 550   fSecRejectRatio(.4),
 551   fROC(AliTPCROC::Instance()),
 552   fMapping(NULL),
 553   fParam(new  AliTPCParam),
 554   fPedestalTPC(0x0),
 555   fPadNoiseTPC(0x0),
 556   fPedestalROC(0x0),
 557   fPadNoiseROC(0x0),
 558   fCalRocArrayT0(72),
 559   fCalRocArrayT0Err(72),
 560   fCalRocArrayQ(72),
 561   fCalRocArrayRMS(72),
 562   fCalRocArrayOutliers(72),
 563   fHistoQArray(72),
 564   fHistoT0Array(72),
 565   fHistoRMSArray(72),
 566   fMeanT0rms(0),
 567   fMeanQrms(0),
 568   fMeanRMSrms(0),
 569   fHistoTmean(72),
 570   fParamArrayEventPol1(72),
 571   fParamArrayEventPol2(72),
 572   fTMeanArrayEvent(72),
 573   fQMeanArrayEvent(72),
 574   fVEventTime(10),
 575   fVEventNumber(10),
 576   fNevents(0),
 577   fTimeStamp(0),
 578   fEventId(-1),
 579   fRunNumber(-1),
 580   fOldRunNumber(-1),
 581   fPadTimesArrayEvent(72),
 582   fPadQArrayEvent(72),
 583   fPadRMSArrayEvent(72),
 584   fPadPedestalArrayEvent(72),
 585   fCurrentChannel(-1),
 586   fCurrentSector(-1),
 587   fCurrentRow(-1),
 588   fMaxPadSignal(-1),
 589   fMaxTimeBin(-1),
 590   fPadSignal(1024),
 591   fPadPedestal(0),
 592   fPadNoise(0),
 593   fVTime0Offset(72),
 594   fVTime0OffsetCounter(72),
 595   fVMeanQ(72),
 596   fVMeanQCounter(72),
 597   fCurrentCETimeRef(0),
 598   //  fEvent(-1),
 599   fDebugStreamer(0x0),
 600   fDebugLevel(0)
 601 {
 602   //
 603   // constructor which uses a tmap as input to set some specific parameters
 604   //
 605   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 606   if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 607   if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
 608   if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
 609   if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
 610   if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
 611   if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
 612   if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
 613   if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
 614   if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
 615   if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
 616   if (config->GetValue("PeakMinus")) fPeakMinus = ((TObjString*)config->GetValue("PeakMinus"))->GetString().Atoi();
 617   if (config->GetValue("PeakPlus")) fPeakPlus = ((TObjString*)config->GetValue("PeakPlus"))->GetString().Atoi();
 618   if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
 619   if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
 620   if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
 621   if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
 622
 623   fParam->Update();
 624 }
 625
 626 //_____________________________________________________________________
 627 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
 628 {
 629   //
 630   // assignment operator
 631   //
 632   if (&source == this) return *this;
 633   new (this) AliTPCCalibCE(source);
 634
 635   return *this;
 636 }
 637 //_____________________________________________________________________
 638 AliTPCCalibCE::~AliTPCCalibCE()
 639 {
 640     //
 641     // destructor
 642     //
 643
 644     fCalRocArrayT0.Delete();
 645     fCalRocArrayT0Err.Delete();
 646     fCalRocArrayQ.Delete();
 647     fCalRocArrayRMS.Delete();
 648     fCalRocArrayOutliers.Delete();
 649
 650     fHistoQArray.Delete();
 651     fHistoT0Array.Delete();
 652     fHistoRMSArray.Delete();
 653
 654     fHistoTmean.Delete();
 655
 656     fParamArrayEventPol1.Delete();
 657     fParamArrayEventPol2.Delete();
 658     fTMeanArrayEvent.Delete();
 659     fQMeanArrayEvent.Delete();
 660
 661     fPadTimesArrayEvent.Delete();
 662     fPadQArrayEvent.Delete();
 663     fPadRMSArrayEvent.Delete();
 664     fPadPedestalArrayEvent.Delete();
 665
 666     if ( fDebugStreamer) delete fDebugStreamer;
 667 //    if ( fHTime0 ) delete fHTime0;
 668 //    delete fROC;
 669     delete fParam;
 670 }
 671 //_____________________________________________________________________
 672 Int_t AliTPCCalibCE::Update(const Int_t icsector,
 673                                 const Int_t icRow,
 674                                 const Int_t icPad,
 675                                 const Int_t icTimeBin,
 676                                 const Float_t csignal)
 677 {
 678   //
 679   // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 680   // no extra analysis necessary. Assumes knowledge of the signal shape!
 681   // assumes that it is looped over consecutive time bins of one pad
 682   //
 683
 684   //temp
 685
 686   if (icRow<0) return 0;
 687   if (icPad<0) return 0;
 688   if (icTimeBin<0) return 0;
 689   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 690
 691   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 692
 693   //init first pad and sector in this event
 694   if ( fCurrentChannel == -1 ) {
 695     fLastSector=-1;
 696     fCurrentChannel = iChannel;
 697     fCurrentSector  = icsector;
 698     fCurrentRow     = icRow;
 699   }
 700
 701   //process last pad if we change to a new one
 702   if ( iChannel != fCurrentChannel ){
 703     ProcessPad();
 704     fLastSector=fCurrentSector;
 705     fCurrentChannel = iChannel;
 706     fCurrentSector  = icsector;
 707     fCurrentRow     = icRow;
 708   }
 709
 710   //fill signals for current pad
 711   fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
 712   if ( csignal > fMaxPadSignal ){
 713     fMaxPadSignal = csignal;
 714     fMaxTimeBin   = icTimeBin;
 715   }
 716   return 0;
 717 }
 718 //_____________________________________________________________________
 719 void AliTPCCalibCE::FindPedestal(Float_t part)
 720 {
 721   //
 722     // find pedestal and noise for the current pad. Use either database or
 723     // truncated mean with part*100%
 724   //
 725   Bool_t noPedestal = kTRUE;
 726
 727     //use pedestal database if set
 728   if (fPedestalTPC&&fPadNoiseTPC){
 729         //only load new pedestals if the sector has changed
 730     if ( fCurrentSector!=fLastSector ){
 731       fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 732       fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 733     }
 734
 735     if ( fPedestalROC&&fPadNoiseROC ){
 736       fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
 737       fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 738       noPedestal   = kFALSE;
 739     }
 740
 741   }
 742
 743     //if we are not running with pedestal database, or for the current sector there is no information
 744     //available, calculate the pedestal and noise on the fly
 745   if ( noPedestal ) {
 746     fPadPedestal = 0;
 747     fPadNoise    = 0;
 748     if ( fIsZeroSuppressed ) return;
 749     const Int_t kPedMax = 100;  //maximum pedestal value
 750     Float_t  max    =  0;
 751     Float_t  maxPos =  0;
 752     Int_t    median =  -1;
 753     Int_t    count0 =  0;
 754     Int_t    count1 =  0;
 755     //
 756     Float_t padSignal=0;
 757     //
 758     UShort_t histo[kPedMax];
 759     memset(histo,0,kPedMax*sizeof(UShort_t));
 760
 761         //fill pedestal histogram
 762     for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 763       padSignal = fPadSignal.GetMatrixArray()[i];
 764       if (padSignal<=0) continue;
 765       if (padSignal>max && i>10) {
 766         max = padSignal;
 767         maxPos = i;
 768       }
 769       if (padSignal>kPedMax-1) continue;
 770       histo[int(padSignal+0.5)]++;
 771       count0++;
 772     }
 773         //find median
 774     for (Int_t i=1; i<kPedMax; ++i){
 775       if (count1<count0*0.5) median=i;
 776       count1+=histo[i];
 777     }
 778         // truncated mean
 779     //
 780     Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 781     //
 782     for (Int_t idelta=1; idelta<10; ++idelta){
 783       if (median-idelta<=0) continue;
 784       if (median+idelta>kPedMax) continue;
 785       if (count<part*count1){
 786         count+=histo[median-idelta];
 787         mean +=histo[median-idelta]*(median-idelta);
 788         rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 789         count+=histo[median+idelta];
 790         mean +=histo[median+idelta]*(median+idelta);
 791         rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 792       }
 793     }
 794     if ( count > 0 ) {
 795       mean/=count;
 796       rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 797       fPadPedestal = mean;
 798       fPadNoise    = rms;
 799     }
 800   }
 801 }
 802 //_____________________________________________________________________
 803 void AliTPCCalibCE::UpdateCETimeRef()
 804 {
 805   // Find the time reference of the last valid CE signal in sector
 806   // for irocs of the A-Side the reference of the corresponging OROC is returned
 807   // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
 808   if ( fLastSector == fCurrentSector ) return;
 809   Int_t sector=fCurrentSector;
 810   if ( sector < 18 ) sector+=36;
 811   fCurrentCETimeRef=0;
 812   TVectorF *vtRef = GetTMeanEvents(sector);
 813   if ( !vtRef ) return;
 814   Int_t vtRefSize= vtRef->GetNrows();
 815   if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
 816   else vtRefSize=fNevents;
 817   while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
 818   fCurrentCETimeRef=(*vtRef)[vtRefSize];
 819   AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
 820 }
 821 //_____________________________________________________________________
 822 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
 823 {
 824   //
 825     //  Find position, signal width and height of the CE signal (last signal)
 826     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 827     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 828   //
 829
 830   Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 831   Int_t   cemaxpos       = 0;
 832   Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise;  // threshold for the signal sum
 833   const Int_t    kCemin  = 4;             // range for the analysis of the ce signal +- channels from the peak
 834   const Int_t    kCemax  = 7;
 835
 836   Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
 837
 838     // find maximum closest to the sector mean from the last event
 839   for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
 840         // get sector mean of last event
 841     Float_t tmean = fCurrentCETimeRef;
 842     if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
 843       minDist  = tmean-maxima[imax];
 844       cemaxpos = (Int_t)maxima[imax];
 845     }
 846   }
 847
 848   if (cemaxpos!=0){
 849     ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 850     for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
 851       if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
 852         Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 853         if (signal>0) {
 854           ceTime+=signal*(i+0.5);
 855           ceRMS +=signal*(i+0.5)*(i+0.5);
 856           ceQsum+=signal;
 857         }
 858       }
 859     }
 860   }
 861   if (ceQmax&&ceQsum>ceSumThreshold) {
 862     ceTime/=ceQsum;
 863     ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 864     fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 865     fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 866
 867   //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
 868   //                                the pick-up signal should scale with the pad area. In addition
 869   //                                the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
 870   //                                ratio 2/3. The pad area we express in cm2. We normalise the signal
 871   //                                to the OROC signal (factor 2/3 for the IROCs).
 872     Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
 873     if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
 874
 875     ceQsum/=norm;
 876     fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
 877     fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
 878   } else {
 879     ceQmax=0;
 880     ceTime=0;
 881     ceRMS =0;
 882     ceQsum=0;
 883   }
 884   param[0] = ceQmax;
 885   param[1] = ceTime;
 886   param[2] = ceRMS;
 887   qSum     = ceQsum;
 888 }
 889 //_____________________________________________________________________
 890 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
 891 {
 892   //
 893     // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
 894     // and 'tplus' timebins after 'pos'
 895   //
 896   if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
 897   for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
 898     if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
 899   for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
 900     if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
 901       ++iTime2;
 902     if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
 903   }
 904   return kTRUE;
 905 }
 906 //_____________________________________________________________________
 907 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
 908 {
 909   //
 910     // Find local maxima on the pad signal and Histogram them
 911   //
 912   Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal
 913   Int_t   count       = 0;
 914 //    Int_t   tminus      = 2;
 915 //    Int_t   tplus       = 3;
 916   for (Int_t i=fLastTimeBin-fPeakPlus-1; i>=fFirstTimeBin+fPeakMinus; --i){
 917     if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,fPeakMinus,fPeakPlus) ){
 918       if (count<maxima.GetNrows()){
 919         maxima.GetMatrixArray()[count++]=i;
 920         GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
 921       }
 922     }
 923   }
 924 }
 925 //_____________________________________________________________________
 926 void AliTPCCalibCE::ProcessPad()
 927 {
 928     //
 929     //  Process data of current pad
 930     //
 931     FindPedestal();
 932
 933     TVectorF maxima(15);     // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
 934                              // + central electrode and possibly post peaks from the CE signal
 935                              // however if we are on a high noise pad a lot more peaks due to the noise might occur
 936     FindLocalMaxima(maxima);
 937     if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
 938
 939     UpdateCETimeRef();                       // update the time refenrence for the current sector
 940     if ( fCurrentCETimeRef==0 ) return;      //return if we don't have time 0 info, eg if only one side has laser
 941     TVectorD param(3);
 942     Float_t  qSum;
 943     FindCESignal(param, qSum, maxima);
 944
 945     Double_t meanT  = param[1];
 946     Double_t sigmaT = param[2];
 947
 948     //Fill Event T0 counter
 949     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 950
 951     //Fill Q histogram
 952     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
 953
 954     //Fill RMS histogram
 955     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 956
 957
 958     //Fill debugging info
 959     if ( fDebugLevel>0 ){
 960         (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
 961         (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
 962         (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
 963     }
 964
 965     ResetPad();
 966 }
 967 //_____________________________________________________________________
 968 void AliTPCCalibCE::EndEvent()
 969 {
 970   //  Process data of current pad
 971   //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
 972   //  before the EndEvent function to set the event timestamp and number!!!
 973   //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
 974   //  function was called
 975
 976   //check if last pad has allready been processed, if not do so
 977   if ( fMaxTimeBin>-1 ) ProcessPad();
 978
 979   AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
 980
 981   TVectorD param(3);
 982   TMatrixD dummy(3,3);
 983 //    TVectorF vMeanTime(72);
 984 //    TVectorF vMeanQ(72);
 985   AliTPCCalROC *calIroc=new AliTPCCalROC(0);
 986   AliTPCCalROC *calOroc=new AliTPCCalROC(36);
 987
 988   //find mean time0 offset for side A and C
 989   //use only orocs due to the better statistics
 990   Double_t time0Side[2];       //time0 for side A:0 and C:1
 991   Double_t time0SideCount[2];  //time0 counter for side A:0 and C:1
 992   time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
 993   for ( Int_t iSec = 36; iSec<72; ++iSec ){
 994     time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
 995     time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
 996   }
 997   if ( time0SideCount[0] >0  )
 998     time0Side[0]/=time0SideCount[0];
 999   if ( time0SideCount[1] >0 )
1000     time0Side[1]/=time0SideCount[1];
1001     // end find time0 offset
1002   AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1003   Int_t nSecMeanT=0;
1004   //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1005   for ( Int_t iSec = 0; iSec<72; ++iSec ){
1006     AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1007     //find median and then calculate the mean around it
1008     TH1S *hMeanT    = GetHistoTmean(iSec); //histogram with local maxima position information
1009     if ( !hMeanT ) continue;
1010     //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1011     if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1012       hMeanT->Reset();
1013       AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1014       continue;
1015     }
1016
1017     Double_t entries = hMeanT->GetEffectiveEntries();
1018     Double_t sum     = 0;
1019     Short_t *arr     = hMeanT->GetArray()+1;
1020     Int_t ibin=0;
1021     for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1022       sum+=arr[ibin];
1023       if ( sum>=(entries/2.) ) break;
1024     }
1025     Int_t delta = 4;
1026     Int_t firstBin = fFirstTimeBin+ibin-delta;
1027     Int_t lastBin  = fFirstTimeBin+ibin+delta;
1028     if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1029     if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
1030     Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1031
1032         // check boundaries for ebye info of mean time
1033     TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1034     Int_t vSize=vMeanTime->GetNrows();
1035     if ( vSize < fNevents+1 ){
1036       vMeanTime->ResizeTo(vSize+100);
1037     }
1038
1039     vMeanTime->GetMatrixArray()[fNevents]=median;
1040     nSecMeanT++;
1041     // end find median
1042
1043     TVectorF *vTimes = GetPadTimesEvent(iSec);
1044     if ( !vTimes ) continue;                     //continue if no time information for this sector is available
1045
1046     AliTPCCalROC calIrocOutliers(0);
1047     AliTPCCalROC calOrocOutliers(36);
1048
1049     // calculate mean Q of the sector
1050     TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1051     vSize=vMeanQ->GetNrows();
1052     if ( vSize < fNevents+1 ){
1053       vMeanQ->ResizeTo(vSize+100);
1054     }
1055     Float_t meanQ = 0;
1056     if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1057     vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1058
1059     for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1060       Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
1061
1062             //set values for temporary roc calibration class
1063       if ( iSec < 36 ) {
1064         calIroc->SetValue(iChannel, time);
1065         if ( time == 0 ) calIrocOutliers.SetValue(iChannel,1);
1066
1067       } else {
1068         calOroc->SetValue(iChannel, time);
1069         if ( time == 0 ) calOrocOutliers.SetValue(iChannel,1);
1070       }
1071
1072       if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1073         GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
1074
1075
1076
1077             //-------------------------------  Debug start  ------------------------------
1078       if ( fDebugLevel>0 ){
1079         if ( !fDebugStreamer ) {
1080                         //debug stream
1081           TDirectory *backup = gDirectory;
1082           fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1083           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1084         }
1085
1086         Int_t row=0;
1087         Int_t pad=0;
1088         Int_t padc=0;
1089
1090         Float_t q   = (*GetPadQEvent(iSec))[iChannel];
1091         Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1092
1093         UInt_t channel=iChannel;
1094         Int_t sector=iSec;
1095
1096         while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1097         pad = channel-fROC->GetRowIndexes(sector)[row];
1098         padc = pad-(fROC->GetNPads(sector,row)/2);
1099
1100 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1101 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
1102 //                                  fLastTimeBin-fFirstTimeBin,
1103 //                                  fFirstTimeBin,fLastTimeBin);
1104 //              h1->SetDirectory(0);
1105         //
1106 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1107 //                  h1->Fill(i,fPadSignal(i));
1108
1109         Double_t t0Sec = 0;
1110         if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1111           t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1112         Double_t t0Side = time0Side[(iSec/18)%2];
1113         (*fDebugStreamer) << "DataPad" <<
1114             "Event=" << fNevents <<
1115             "RunNumber=" << fRunNumber <<
1116             "TimeStamp="   << fTimeStamp <<
1117             "Sector="<< sector <<
1118             "Row="   << row<<
1119             "Pad="   << pad <<
1120             "PadC="  << padc <<
1121             "PadSec="<< channel <<
1122             "Time0Sec="  << t0Sec <<
1123             "Time0Side=" << t0Side <<
1124             "Time="  << time <<
1125             "RMS="   << rms <<
1126             "Sum="   << q <<
1127             "MeanQ=" << meanQ <<
1128                     //              "hist.=" << h1 <<
1129             "\n";
1130
1131                 //              delete h1;
1132
1133       }
1134             //-----------------------------  Debug end  ------------------------------
1135     }// end channel loop
1136
1137     TVectorD paramPol1(3);
1138     TVectorD paramPol2(6);
1139     TMatrixD matPol1(3,3);
1140     TMatrixD matPol2(6,6);
1141     Float_t  chi2Pol1=0;
1142     Float_t  chi2Pol2=0;
1143
1144     if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1145       if ( iSec < 36 ){
1146         calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1147         calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1148       } else {
1149         calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1150         calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1151       }
1152
1153       GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1154       GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1155     }
1156
1157         //-------------------------------  Debug start  ------------------------------
1158     if ( fDebugLevel>0 ){
1159       if ( !fDebugStreamer ) {
1160                 //debug stream
1161         TDirectory *backup = gDirectory;
1162         fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1163         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1164       }
1165       (*fDebugStreamer) << "DataRoc" <<
1166 //              "Event=" << fEvent <<
1167           "RunNumber=" << fRunNumber <<
1168           "TimeStamp="   << fTimeStamp <<
1169           "Sector="<< iSec <<
1170           "hMeanT.=" << hMeanT <<
1171           "median=" << median <<
1172           "paramPol1.=" << &paramPol1 <<
1173           "paramPol2.=" << &paramPol2 <<
1174           "matPol1.="   << &matPol1 <<
1175           "matPol2.="   << &matPol2 <<
1176           "chi2Pol1="   << chi2Pol1 <<
1177           "chi2Pol2="   << chi2Pol2 <<
1178           "\n";
1179     }
1180         //-------------------------------  Debug end  ------------------------------
1181     hMeanT->Reset();
1182   }// end sector loop
1183     //return if no sector has a valid mean time
1184   if ( nSecMeanT == 0 ) return;
1185
1186
1187 //    fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1188 //    fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1189   if ( fVEventTime.GetNrows() < fNevents+1 ) {
1190     fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1191     fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1192   }
1193   fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1194   fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1195
1196   fNevents++;
1197   fOldRunNumber = fRunNumber;
1198
1199   delete calIroc;
1200   delete calOroc;
1201   AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1202 }
1203 //_____________________________________________________________________
1204 Bool_t AliTPCCalibCE::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
1205 {
1206   //
1207   // Event Processing loop - AliTPCRawStreamFast
1208   //
1209   ResetEvent();
1210   Bool_t withInput = kFALSE;
1211   while ( rawStreamFast->NextDDL() ){
1212     while ( rawStreamFast->NextChannel() ){
1213       Int_t isector  = rawStreamFast->GetSector();                       //  current sector
1214       Int_t iRow     = rawStreamFast->GetRow();                          //  current row
1215       Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
1216
1217       while ( rawStreamFast->NextBunch() ){
1218         Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
1219         Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
1220         for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
1221           Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
1222           Update(isector,iRow,iPad,iTimeBin+1,signal);
1223           withInput = kTRUE;
1224         }
1225       }
1226     }
1227   }
1228   if (withInput){
1229     EndEvent();
1230   }
1231   return withInput;
1232 }
1233 //_____________________________________________________________________
1234 Bool_t AliTPCCalibCE::ProcessEventFast(AliRawReader *rawReader)
1235 {
1236   //
1237   //  Event processing loop using the fast raw stream algorithm- AliRawReader
1238   //
1239
1240   //printf("ProcessEventFast - raw reader\n");
1241
1242   AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1243   if (eventHeader){
1244       fTimeStamp   = eventHeader->Get("Timestamp");
1245       fRunNumber = eventHeader->Get("RunNb");
1246   }
1247   fEventId = *rawReader->GetEventId();
1248
1249   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
1250   Bool_t res=ProcessEventFast(rawStreamFast);
1251   delete rawStreamFast;
1252   return res;
1253
1254 }
1255 //_____________________________________________________________________
1256 Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream)
1257 {
1258   //
1259   // Event Processing loop - AliTPCRawStream
1260   // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
1261   //
1262
1263   ResetEvent();
1264
1265   Bool_t withInput = kFALSE;
1266
1267   while (rawStream->Next()) {
1268
1269       Int_t isector  = rawStream->GetSector();                       //  current sector
1270       Int_t iRow     = rawStream->GetRow();                          //  current row
1271       Int_t iPad     = rawStream->GetPad();                          //  current pad
1272       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
1273       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
1274
1275       Update(isector,iRow,iPad,iTimeBin,signal);
1276       withInput = kTRUE;
1277   }
1278
1279   if (withInput){
1280       EndEvent();
1281   }
1282
1283   return withInput;
1284 }
1285 //_____________________________________________________________________
1286 Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
1287 {
1288   //
1289   //  Event processing loop - AliRawReader
1290   //
1291
1292
1293   AliTPCRawStream rawStream(rawReader,(AliAltroMapping**)fMapping);
1294     AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1295     if (eventHeader){
1296         fTimeStamp   = eventHeader->Get("Timestamp");
1297         fRunNumber = eventHeader->Get("RunNb");
1298     }
1299     fEventId = *rawReader->GetEventId();
1300
1301
1302     rawReader->Select("TPC");
1303
1304   return ProcessEvent(&rawStream);
1305 }
1306 //_____________________________________________________________________
1307 Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
1308 {
1309   //
1310   //  Event processing loop - date event
1311   //
1312     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1313     Bool_t result=ProcessEvent(rawReader);
1314     delete rawReader;
1315     return result;
1316
1317 }
1318 //_____________________________________________________________________
1319 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1320                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
1321                                   Char_t *type, Bool_t force)
1322 {
1323     //
1324     // return pointer to TH2S histogram of 'type'
1325     // if force is true create a new histogram if it doesn't exist allready
1326     //
1327     if ( !force || arr->UncheckedAt(sector) )
1328         return (TH2S*)arr->UncheckedAt(sector);
1329
1330     // if we are forced and histogram doesn't exist yet create it
1331     Char_t name[255], title[255];
1332
1333     sprintf(name,"hCalib%s%.2d",type,sector);
1334     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1335
1336     // new histogram with Q calib information. One value for each pad!
1337     TH2S* hist = new TH2S(name,title,
1338                           nbinsY, ymin, ymax,
1339                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1340     hist->SetDirectory(0);
1341     arr->AddAt(hist,sector);
1342     return hist;
1343 }
1344 //_____________________________________________________________________
1345 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1346 {
1347     //
1348     // return pointer to T0 histogram
1349     // if force is true create a new histogram if it doesn't exist allready
1350     //
1351     TObjArray *arr = &fHistoT0Array;
1352     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1353 }
1354 //_____________________________________________________________________
1355 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1356 {
1357     //
1358     // return pointer to Q histogram
1359     // if force is true create a new histogram if it doesn't exist allready
1360     //
1361     TObjArray *arr = &fHistoQArray;
1362     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1363 }
1364 //_____________________________________________________________________
1365 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1366 {
1367     //
1368     // return pointer to Q histogram
1369     // if force is true create a new histogram if it doesn't exist allready
1370     //
1371     TObjArray *arr = &fHistoRMSArray;
1372     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1373 }
1374 //_____________________________________________________________________
1375 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1376                               Char_t *type, Bool_t force)
1377 {
1378     //
1379     // return pointer to TH1S histogram
1380     // if force is true create a new histogram if it doesn't exist allready
1381     //
1382     if ( !force || arr->UncheckedAt(sector) )
1383         return (TH1S*)arr->UncheckedAt(sector);
1384
1385     // if we are forced and histogram doesn't yes exist create it
1386     Char_t name[255], title[255];
1387
1388     sprintf(name,"hCalib%s%.2d",type,sector);
1389     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1390
1391     // new histogram with calib information. One value for each pad!
1392     TH1S* hist = new TH1S(name,title,
1393                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1394     hist->SetDirectory(0);
1395     arr->AddAt(hist,sector);
1396     return hist;
1397 }
1398 //_____________________________________________________________________
1399 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1400 {
1401     //
1402     // return pointer to Q histogram
1403     // if force is true create a new histogram if it doesn't exist allready
1404     //
1405     TObjArray *arr = &fHistoTmean;
1406     return GetHisto(sector, arr, "LastTmean", force);
1407 }
1408 //_____________________________________________________________________
1409 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1410 {
1411     //
1412     // return pointer to Pad Info from 'arr' for the current event and sector
1413     // if force is true create it if it doesn't exist allready
1414     //
1415     if ( !force || arr->UncheckedAt(sector) )
1416         return (TVectorF*)arr->UncheckedAt(sector);
1417
1418     TVectorF *vect = new TVectorF(size);
1419     arr->AddAt(vect,sector);
1420     return vect;
1421 }
1422 //_____________________________________________________________________
1423 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1424 {
1425     //
1426     // return pointer to Pad Times Array for the current event and sector
1427     // if force is true create it if it doesn't exist allready
1428     //
1429     TObjArray *arr = &fPadTimesArrayEvent;
1430     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1431 }
1432 //_____________________________________________________________________
1433 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1434 {
1435     //
1436     // return pointer to Pad Q Array for the current event and sector
1437     // if force is true create it if it doesn't exist allready
1438     // for debugging purposes only
1439     //
1440
1441     TObjArray *arr = &fPadQArrayEvent;
1442     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1443 }
1444 //_____________________________________________________________________
1445 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1446 {
1447     //
1448     // return pointer to Pad RMS Array for the current event and sector
1449     // if force is true create it if it doesn't exist allready
1450     // for debugging purposes only
1451     //
1452     TObjArray *arr = &fPadRMSArrayEvent;
1453     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1454 }
1455 //_____________________________________________________________________
1456 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1457 {
1458     //
1459     // return pointer to Pad RMS Array for the current event and sector
1460     // if force is true create it if it doesn't exist allready
1461     // for debugging purposes only
1462     //
1463     TObjArray *arr = &fPadPedestalArrayEvent;
1464     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1465 }
1466 //_____________________________________________________________________
1467 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1468 {
1469     //
1470     // return pointer to the EbyE info of the mean arrival time for 'sector'
1471     // if force is true create it if it doesn't exist allready
1472     //
1473     TObjArray *arr = &fTMeanArrayEvent;
1474     return GetVectSector(sector,arr,100,force);
1475 }
1476 //_____________________________________________________________________
1477 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1478 {
1479     //
1480     // return pointer to the EbyE info of the mean arrival time for 'sector'
1481     // if force is true create it if it doesn't exist allready
1482     //
1483     TObjArray *arr = &fQMeanArrayEvent;
1484     return GetVectSector(sector,arr,100,force);
1485 }
1486 //_____________________________________________________________________
1487 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1488 {
1489     //
1490     // return pointer to ROC Calibration
1491     // if force is true create a new histogram if it doesn't exist allready
1492     //
1493     if ( !force || arr->UncheckedAt(sector) )
1494         return (AliTPCCalROC*)arr->UncheckedAt(sector);
1495
1496     // if we are forced and histogram doesn't yes exist create it
1497
1498     // new AliTPCCalROC for T0 information. One value for each pad!
1499     AliTPCCalROC *croc = new AliTPCCalROC(sector);
1500     arr->AddAt(croc,sector);
1501     return croc;
1502 }
1503 //_____________________________________________________________________
1504 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1505 {
1506     //
1507     // return pointer to Time 0 ROC Calibration
1508     // if force is true create a new histogram if it doesn't exist allready
1509     //
1510     TObjArray *arr = &fCalRocArrayT0;
1511     return GetCalRoc(sector, arr, force);
1512 }
1513 //_____________________________________________________________________
1514 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1515 {
1516     //
1517     // return pointer to the error of Time 0 ROC Calibration
1518     // if force is true create a new histogram if it doesn't exist allready
1519     //
1520     TObjArray *arr = &fCalRocArrayT0Err;
1521     return GetCalRoc(sector, arr, force);
1522 }
1523 //_____________________________________________________________________
1524 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1525 {
1526     //
1527     // return pointer to T0 ROC Calibration
1528     // if force is true create a new histogram if it doesn't exist allready
1529     //
1530     TObjArray *arr = &fCalRocArrayQ;
1531     return GetCalRoc(sector, arr, force);
1532 }
1533 //_____________________________________________________________________
1534 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1535 {
1536     //
1537     // return pointer to signal width ROC Calibration
1538     // if force is true create a new histogram if it doesn't exist allready
1539     //
1540     TObjArray *arr = &fCalRocArrayRMS;
1541     return GetCalRoc(sector, arr, force);
1542 }
1543 //_____________________________________________________________________
1544 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1545 {
1546     //
1547     // return pointer to Outliers
1548     // if force is true create a new histogram if it doesn't exist allready
1549     //
1550     TObjArray *arr = &fCalRocArrayOutliers;
1551     return GetCalRoc(sector, arr, force);
1552 }
1553 //_____________________________________________________________________
1554 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1555 {
1556     //
1557     // return pointer to TObjArray of fit parameters
1558     // if force is true create a new histogram if it doesn't exist allready
1559     //
1560     if ( !force || arr->UncheckedAt(sector) )
1561         return (TObjArray*)arr->UncheckedAt(sector);
1562
1563     // if we are forced and array doesn't yes exist create it
1564
1565     // new TObjArray for parameters
1566     TObjArray *newArr = new TObjArray;
1567     arr->AddAt(newArr,sector);
1568     return newArr;
1569 }
1570 //_____________________________________________________________________
1571 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1572 {
1573     //
1574     // return pointer to TObjArray of fit parameters from plane fit
1575     // if force is true create a new histogram if it doesn't exist allready
1576     //
1577     TObjArray *arr = &fParamArrayEventPol1;
1578     return GetParamArray(sector, arr, force);
1579 }
1580 //_____________________________________________________________________
1581 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1582 {
1583     //
1584     // return pointer to TObjArray of fit parameters from parabola fit
1585     // if force is true create a new histogram if it doesn't exist allready
1586     //
1587     TObjArray *arr = &fParamArrayEventPol2;
1588     return GetParamArray(sector, arr, force);
1589 }
1590 //_____________________________________________________________________
1591 void AliTPCCalibCE::ResetEvent()
1592 {
1593     //
1594     //  Reset global counters  -- Should be called before each event is processed
1595     //
1596     fLastSector=-1;
1597     fCurrentSector=-1;
1598     fCurrentRow=-1;
1599     fCurrentChannel=-1;
1600
1601     ResetPad();
1602
1603     fPadTimesArrayEvent.Delete();
1604     fPadQArrayEvent.Delete();
1605     fPadRMSArrayEvent.Delete();
1606     fPadPedestalArrayEvent.Delete();
1607
1608     for ( Int_t i=0; i<72; ++i ){
1609         fVTime0Offset.GetMatrixArray()[i]=0;
1610         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1611         fVMeanQ.GetMatrixArray()[i]=0;
1612         fVMeanQCounter.GetMatrixArray()[i]=0;
1613     }
1614 }
1615 //_____________________________________________________________________
1616 void AliTPCCalibCE::ResetPad()
1617 {
1618     //
1619     //  Reset pad infos -- Should be called after a pad has been processed
1620     //
1621     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1622         fPadSignal.GetMatrixArray()[i] = 0;
1623     fMaxTimeBin   = -1;
1624     fMaxPadSignal = -1;
1625     fPadPedestal  = -1;
1626     fPadNoise     = -1;
1627 }
1628 //_____________________________________________________________________
1629 void AliTPCCalibCE::Merge(AliTPCCalibCE *ce)
1630 {
1631     //
1632     //  Merge ce to the current AliTPCCalibCE
1633     //
1634
1635     //merge histograms
1636     for (Int_t iSec=0; iSec<72; ++iSec){
1637         TH2S *hRefQmerge   = ce->GetHistoQ(iSec);
1638         TH2S *hRefT0merge  = ce->GetHistoT0(iSec);
1639         TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1640
1641
1642         if ( hRefQmerge ){
1643             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1644             TH2S *hRefQ   = GetHistoQ(iSec);
1645             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1646             else {
1647                 TH2S *hist = new TH2S(*hRefQmerge);
1648                 hist->SetDirectory(0);
1649                 fHistoQArray.AddAt(hist, iSec);
1650             }
1651             hRefQmerge->SetDirectory(dir);
1652         }
1653         if ( hRefT0merge ){
1654             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1655             TH2S *hRefT0  = GetHistoT0(iSec);
1656             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1657             else {
1658                 TH2S *hist = new TH2S(*hRefT0merge);
1659                 hist->SetDirectory(0);
1660                 fHistoT0Array.AddAt(hist, iSec);
1661             }
1662             hRefT0merge->SetDirectory(dir);
1663         }
1664         if ( hRefRMSmerge ){
1665             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1666             TH2S *hRefRMS = GetHistoRMS(iSec);
1667             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1668             else {
1669                 TH2S *hist = new TH2S(*hRefRMSmerge);
1670                 hist->SetDirectory(0);
1671                 fHistoRMSArray.AddAt(hist, iSec);
1672             }
1673             hRefRMSmerge->SetDirectory(dir);
1674         }
1675
1676     }
1677
1678     // merge time information
1679
1680
1681     Int_t nCEevents = ce->GetNeventsProcessed();
1682     for (Int_t iSec=0; iSec<72; ++iSec){
1683         TObjArray *arrPol1CE  = ce->GetParamArrayPol1(iSec);
1684         TObjArray *arrPol2CE  = ce->GetParamArrayPol2(iSec);
1685         TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1686         TVectorF *vMeanQCE    = ce->GetQMeanEvents(iSec);
1687
1688         TObjArray *arrPol1  = 0x0;
1689         TObjArray *arrPol2  = 0x0;
1690         TVectorF *vMeanTime = 0x0;
1691         TVectorF *vMeanQ    = 0x0;
1692
1693         //resize arrays
1694         if ( arrPol1CE && arrPol2CE ){
1695             arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1696             arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1697             arrPol1->Expand(fNevents+nCEevents);
1698             arrPol2->Expand(fNevents+nCEevents);
1699         }
1700         if ( vMeanTimeCE && vMeanQCE ){
1701             vMeanTime = GetTMeanEvents(iSec,kTRUE);
1702             vMeanQ    = GetQMeanEvents(iSec,kTRUE);
1703             vMeanTime->ResizeTo(fNevents+nCEevents);
1704             vMeanQ->ResizeTo(fNevents+nCEevents);
1705         }
1706
1707
1708         for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1709             if ( arrPol1CE && arrPol2CE ){
1710                 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1711                 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1712                 if ( paramPol1 && paramPol2 ){
1713                     GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1714                     GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1715                 }
1716             }
1717             if ( vMeanTimeCE && vMeanQCE ){
1718                 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1719                 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1720             }
1721         }
1722     }
1723
1724
1725
1726     TVectorD*  eventTimes  = ce->GetEventTimes();
1727     TVectorD*  eventIds  = ce->GetEventIds();
1728     fVEventTime.ResizeTo(fNevents+nCEevents);
1729     fVEventNumber.ResizeTo(fNevents+nCEevents);
1730
1731     for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1732         Double_t evTime     = eventTimes->GetMatrixArray()[iEvent];
1733         Double_t evId       = eventIds->GetMatrixArray()[iEvent];
1734         fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1735         fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1736     }
1737     fNevents+=nCEevents; //increase event counter
1738
1739 }
1740 //_____________________________________________________________________
1741 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1742 {
1743   //
1744   // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1745   // xVariable:    0-event time, 1-event id, 2-internal event counter
1746   // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1747   // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1748   //                               0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1749   //                               not used for mean time and mean Q )
1750   // for an example see class description at the beginning
1751   //
1752
1753   Double_t *x = new Double_t[fNevents];
1754   Double_t *y = new Double_t[fNevents];
1755
1756   TVectorD *xVar = 0x0;
1757   TObjArray *aType = 0x0;
1758   Int_t npoints=0;
1759
1760     // sanity checks
1761   if ( !GetHistoT0(sector) )            return 0x0; //Sector has not been filled
1762   if ( (sector<0) || (sector>71) )      return 0x0;
1763   if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1764   if ( (fitType<0) || (fitType>3) )     return 0x0;
1765   if ( !GetTMeanEvents(sector) )        return 0x0; //no mean time information available
1766
1767   if ( fitType==0 ){
1768     if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1769     aType = &fParamArrayEventPol1;
1770     if ( aType->At(sector)==0x0 ) return 0x0;
1771   }
1772   else if ( fitType==1 ){
1773     if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1774     aType = &fParamArrayEventPol2;
1775     if ( aType->At(sector)==0x0 ) return 0x0;
1776   }
1777
1778
1779   if ( xVariable == 0 ) xVar = &fVEventTime;
1780   if ( xVariable == 1 ) xVar = &fVEventNumber;
1781   if ( xVariable == 2 ) {
1782     xVar = new TVectorD(fNevents);
1783     for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1784   }
1785
1786   for (Int_t ievent =0; ievent<fNevents; ++ievent){
1787     if ( fitType<2 ){
1788       TObjArray *events = (TObjArray*)(aType->At(sector));
1789       if ( events->GetSize()<=ievent ) break;
1790       TVectorD *v = (TVectorD*)(events->At(ievent));
1791       if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1792     } else if (fitType == 2) {
1793       Double_t xValue=(*xVar)[ievent];
1794       Double_t yValue=(*GetTMeanEvents(sector))[ievent];
1795       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1796     }else if (fitType == 3) {
1797       Double_t xValue=(*xVar)[ievent];
1798       Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1799       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1800     }
1801   }
1802
1803   TGraph *gr = new TGraph(npoints);
1804     //sort xVariable increasing
1805   Int_t    *sortIndex = new Int_t[npoints];
1806   TMath::Sort(npoints,x,sortIndex);
1807   for (Int_t i=0;i<npoints;++i){
1808     gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1809   }
1810
1811
1812   if ( xVariable == 2 ) delete xVar;
1813   delete x;
1814   delete y;
1815   delete sortIndex;
1816   return gr;
1817 }
1818 //_____________________________________________________________________
1819 void AliTPCCalibCE::Analyse()
1820 {
1821     //
1822     //  Calculate calibration constants
1823     //
1824
1825     TVectorD paramQ(3);
1826     TVectorD paramT0(3);
1827     TVectorD paramRMS(3);
1828     TMatrixD dummy(3,3);
1829
1830     Float_t channelCounter=0;
1831     fMeanT0rms=0;
1832     fMeanQrms=0;
1833     fMeanRMSrms=0;
1834
1835     for (Int_t iSec=0; iSec<72; ++iSec){
1836         TH2S *hT0 = GetHistoT0(iSec);
1837         if (!hT0 ) continue;
1838
1839         AliTPCCalROC *rocQ     = GetCalRocQ  (iSec,kTRUE);
1840         AliTPCCalROC *rocT0    = GetCalRocT0 (iSec,kTRUE);
1841         AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1842         AliTPCCalROC *rocRMS   = GetCalRocRMS(iSec,kTRUE);
1843         AliTPCCalROC *rocOut   = GetCalRocOutliers(iSec,kTRUE);
1844
1845         TH2S *hQ   = GetHistoQ(iSec);
1846         TH2S *hRMS = GetHistoRMS(iSec);
1847
1848         Short_t *arrayhQ   = hQ->GetArray();
1849         Short_t *arrayhT0  = hT0->GetArray();
1850         Short_t *arrayhRMS = hRMS->GetArray();
1851
1852         UInt_t nChannels = fROC->GetNChannels(iSec);
1853
1854         //debug
1855         Int_t row=0;
1856         Int_t pad=0;
1857         Int_t padc=0;
1858         //! debug
1859
1860         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1861
1862
1863             Float_t cogTime0 = -1000;
1864             Float_t cogQ     = -1000;
1865             Float_t cogRMS   = -1000;
1866             Float_t cogOut   = 0;
1867             Float_t rms      = 0;
1868             Float_t rmsT0    = 0;
1869
1870
1871             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1872             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1873             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1874
1875             cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
1876             fMeanQrms+=rms;
1877             cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
1878             fMeanT0rms+=rmsT0;
1879             cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
1880             fMeanRMSrms+=rms;
1881             channelCounter++;
1882
1883             /*
1884              //outlier specifications
1885             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1886                 cogOut = 1;
1887                 cogTime0 = 0;
1888                 cogQ     = 0;
1889                 cogRMS   = 0;
1890             }
1891 */
1892             rocQ->SetValue(iChannel, cogQ*cogQ);
1893             rocT0->SetValue(iChannel, cogTime0);
1894             rocT0Err->SetValue(iChannel, rmsT0);
1895             rocRMS->SetValue(iChannel, cogRMS);
1896             rocOut->SetValue(iChannel, cogOut);
1897
1898
1899             //debug
1900             if ( fDebugLevel > 0 ){
1901                 if ( !fDebugStreamer ) {
1902                         //debug stream
1903                     TDirectory *backup = gDirectory;
1904                     fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1905                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1906                 }
1907
1908                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1909                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1910                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1911
1912                 (*fDebugStreamer) << "DataEnd" <<
1913                     "Sector="  << iSec      <<
1914                     "Pad="     << pad       <<
1915                     "PadC="    << padc      <<
1916                     "Row="     << row       <<
1917                     "PadSec="  << iChannel   <<
1918                     "Q="       << cogQ      <<
1919                     "T0="      << cogTime0  <<
1920                     "RMS="     << cogRMS    <<
1921                     "\n";
1922             }
1923             //! debug
1924
1925         }
1926
1927     }
1928     if ( channelCounter>0 ){
1929         fMeanT0rms/=channelCounter;
1930         fMeanQrms/=channelCounter;
1931         fMeanRMSrms/=channelCounter;
1932     }
1933     if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
1934 //    delete fDebugStreamer;
1935 //    fDebugStreamer = 0x0;
1936 }
1937 //_____________________________________________________________________
1938 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1939 {
1940     //
1941     //  Write class to file
1942     //
1943
1944     TString sDir(dir);
1945     TString option;
1946
1947     if ( append )
1948         option = "update";
1949     else
1950         option = "recreate";
1951
1952     TDirectory *backup = gDirectory;
1953     TFile f(filename,option.Data());
1954     f.cd();
1955     if ( !sDir.IsNull() ){
1956         f.mkdir(sDir.Data());
1957         f.cd(sDir);
1958     }
1959     this->Write();
1960     f.Close();
1961
1962     if ( backup ) backup->cd();
1963 }