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
 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     fROC(AliTPCROC::Instance()),
 319     fMapping(NULL),
 320     fParam(new AliTPCParam),
 321     fPedestalTPC(0x0),
 322     fPadNoiseTPC(0x0),
 323     fPedestalROC(0x0),
 324     fPadNoiseROC(0x0),
 325     fCalRocArrayT0(72),
 326     fCalRocArrayT0Err(72),
 327     fCalRocArrayQ(72),
 328     fCalRocArrayRMS(72),
 329     fCalRocArrayOutliers(72),
 330     fHistoQArray(72),
 331     fHistoT0Array(72),
 332     fHistoRMSArray(72),
 333     fMeanT0rms(0),
 334     fMeanQrms(0),
 335     fMeanRMSrms(0),
 336     fHistoTmean(72),
 337     fParamArrayEventPol1(72),
 338     fParamArrayEventPol2(72),
 339     fTMeanArrayEvent(72),
 340     fQMeanArrayEvent(72),
 341     fVEventTime(10),
 342     fVEventNumber(10),
 343     fNevents(0),
 344     fTimeStamp(0),
 345     fEventId(-1),
 346     fRunNumber(-1),
 347     fOldRunNumber(-1),
 348     fPadTimesArrayEvent(72),
 349     fPadQArrayEvent(72),
 350     fPadRMSArrayEvent(72),
 351     fPadPedestalArrayEvent(72),
 352     fCurrentChannel(-1),
 353     fCurrentSector(-1),
 354     fCurrentRow(-1),
 355     fMaxPadSignal(-1),
 356     fMaxTimeBin(-1),
 357     fPadSignal(1024),
 358     fPadPedestal(0),
 359     fPadNoise(0),
 360     fVTime0Offset(72),
 361     fVTime0OffsetCounter(72),
 362     fVMeanQ(72),
 363     fVMeanQCounter(72),
 364 //    fEvent(-1),
 365     fDebugStreamer(0x0),
 366     fDebugLevel(0)
 367 {
 368     //
 369     // AliTPCSignal default constructor
 370     //
 371 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
 372     fParam->Update();
 373 }
 374 //_____________________________________________________________________
 375 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
 376     TObject(sig),
 377     fFirstTimeBin(sig.fFirstTimeBin),
 378     fLastTimeBin(sig.fLastTimeBin),
 379     fNbinsT0(sig.fNbinsT0),
 380     fXminT0(sig.fXminT0),
 381     fXmaxT0(sig.fXmaxT0),
 382     fNbinsQ(sig.fNbinsQ),
 383     fXminQ(sig.fXminQ),
 384     fXmaxQ(sig.fXmaxQ),
 385     fNbinsRMS(sig.fNbinsRMS),
 386     fXminRMS(sig.fXminRMS),
 387     fXmaxRMS(sig.fXmaxRMS),
 388     fPeakMinus(sig.fPeakMinus),
 389     fPeakPlus(sig.fPeakPlus),
 390     fNoiseThresholdMax(sig.fNoiseThresholdMax),
 391     fNoiseThresholdSum(sig.fNoiseThresholdSum),
 392     fIsZeroSuppressed(sig.fIsZeroSuppressed),
 393     fLastSector(-1),
 394     fROC(AliTPCROC::Instance()),
 395     fMapping(NULL),
 396     fParam(new AliTPCParam),
 397     fPedestalTPC(0x0),
 398     fPadNoiseTPC(0x0),
 399     fPedestalROC(0x0),
 400     fPadNoiseROC(0x0),
 401     fCalRocArrayT0(72),
 402     fCalRocArrayT0Err(72),
 403     fCalRocArrayQ(72),
 404     fCalRocArrayRMS(72),
 405     fCalRocArrayOutliers(72),
 406     fHistoQArray(72),
 407     fHistoT0Array(72),
 408     fHistoRMSArray(72),
 409     fMeanT0rms(sig.fMeanT0rms),
 410     fMeanQrms(sig.fMeanQrms),
 411     fMeanRMSrms(sig.fMeanRMSrms),
 412     fHistoTmean(72),
 413     fParamArrayEventPol1(72),
 414     fParamArrayEventPol2(72),
 415     fTMeanArrayEvent(72),
 416     fQMeanArrayEvent(72),
 417     fVEventTime(1000),
 418     fVEventNumber(1000),
 419     fNevents(sig.fNevents),
 420     fTimeStamp(0),
 421     fEventId(-1),
 422     fRunNumber(-1),
 423     fOldRunNumber(-1),
 424     fPadTimesArrayEvent(72),
 425     fPadQArrayEvent(72),
 426     fPadRMSArrayEvent(72),
 427     fPadPedestalArrayEvent(72),
 428     fCurrentChannel(-1),
 429     fCurrentSector(-1),
 430     fCurrentRow(-1),
 431     fMaxPadSignal(-1),
 432     fMaxTimeBin(-1),
 433     fPadSignal(1024),
 434     fPadPedestal(0),
 435     fPadNoise(0),
 436     fVTime0Offset(72),
 437     fVTime0OffsetCounter(72),
 438     fVMeanQ(72),
 439     fVMeanQCounter(72),
 440 //    fEvent(-1),
 441     fDebugStreamer(0x0),
 442     fDebugLevel(sig.fDebugLevel)
 443 {
 444     //
 445     // AliTPCSignal copy constructor
 446     //
 447
 448     for (Int_t iSec = 0; iSec < 72; ++iSec){
 449         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 450         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 451         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 452         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 453
 454         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 455         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 456         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 457
 458         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 459         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 460         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 461         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 462
 463         if ( hQ != 0x0 ){
 464 //          TDirectory *dir = hQ->GetDirectory();
 465 //          hQ->SetDirectory(0);
 466             TH2S *hNew = new TH2S(*hQ);
 467             hNew->SetDirectory(0);
 468             fHistoQArray.AddAt(hNew,iSec);
 469 //            hQ->SetDirectory(dir);
 470         }
 471         if ( hT0 != 0x0 ){
 472 //          TDirectory *dir = hT0->GetDirectory();
 473 //          hT0->SetDirectory(0);
 474             TH2S *hNew = new TH2S(*hT0);
 475             hNew->SetDirectory(0);
 476             fHistoT0Array.AddAt(hNew,iSec);
 477 //            hT0->SetDirectory(dir);
 478         }
 479         if ( hRMS != 0x0 ){
 480 //          TDirectory *dir = hRMS->GetDirectory();
 481 //          hRMS->SetDirectory(0);
 482             TH2S *hNew = new TH2S(*hRMS);
 483             hNew->SetDirectory(0);
 484             fHistoRMSArray.AddAt(hNew,iSec);
 485 //            hRMS->SetDirectory(dir);
 486         }
 487     }
 488
 489     //copy fit parameters event by event
 490     TObjArray *arr=0x0;
 491     for (Int_t iSec=0; iSec<72; ++iSec){
 492         arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
 493         if ( arr ){
 494             TObjArray *arrEvents = new TObjArray(arr->GetSize());
 495             fParamArrayEventPol1.AddAt(arrEvents, iSec);
 496             for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 497                 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 498                     arrEvents->AddAt(new TVectorD(*vec),iEvent);
 499         }
 500
 501         arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
 502         if ( arr ){
 503             TObjArray *arrEvents = new TObjArray(arr->GetSize());
 504             fParamArrayEventPol2.AddAt(arrEvents, iSec);
 505             for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 506                 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 507                     arrEvents->AddAt(new TVectorD(*vec),iEvent);
 508         }
 509
 510         TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
 511         TVectorF *vMeanQ    = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
 512         if ( vMeanTime )
 513             fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
 514         if ( vMeanQ )
 515             fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
 516     }
 517
 518
 519     fVEventTime.ResizeTo(sig.fVEventTime);
 520     fVEventNumber.ResizeTo(sig.fVEventNumber);
 521     fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
 522     fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
 523
 524     fParam->Update();
 525 }
 526 //_____________________________________________________________________
 527 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
 528 {
 529   //
 530   // assignment operator
 531   //
 532   if (&source == this) return *this;
 533   new (this) AliTPCCalibCE(source);
 534
 535   return *this;
 536 }
 537 //_____________________________________________________________________
 538 AliTPCCalibCE::~AliTPCCalibCE()
 539 {
 540     //
 541     // destructor
 542     //
 543
 544     fCalRocArrayT0.Delete();
 545     fCalRocArrayT0Err.Delete();
 546     fCalRocArrayQ.Delete();
 547     fCalRocArrayRMS.Delete();
 548     fCalRocArrayOutliers.Delete();
 549
 550     fHistoQArray.Delete();
 551     fHistoT0Array.Delete();
 552     fHistoRMSArray.Delete();
 553
 554     fHistoTmean.Delete();
 555
 556     fParamArrayEventPol1.Delete();
 557     fParamArrayEventPol2.Delete();
 558     fTMeanArrayEvent.Delete();
 559     fQMeanArrayEvent.Delete();
 560
 561     fPadTimesArrayEvent.Delete();
 562     fPadQArrayEvent.Delete();
 563     fPadRMSArrayEvent.Delete();
 564     fPadPedestalArrayEvent.Delete();
 565
 566     if ( fDebugStreamer) delete fDebugStreamer;
 567 //    if ( fHTime0 ) delete fHTime0;
 568 //    delete fROC;
 569     delete fParam;
 570 }
 571 //_____________________________________________________________________
 572 Int_t AliTPCCalibCE::Update(const Int_t icsector,
 573                                 const Int_t icRow,
 574                                 const Int_t icPad,
 575                                 const Int_t icTimeBin,
 576                                 const Float_t csignal)
 577 {
 578     //
 579     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 580     // no extra analysis necessary. Assumes knowledge of the signal shape!
 581     // assumes that it is looped over consecutive time bins of one pad
 582     //
 583
 584     //temp
 585 //    if (icsector<36) return 0;
 586 //    if (icsector%36>17) return 0;
 587
 588
 589   if (icRow<0) return 0;
 590   if (icPad<0) return 0;
 591   if (icTimeBin<0) return 0;
 592     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 593
 594     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 595
 596     //init first pad and sector in this event
 597     if ( fCurrentChannel == -1 ) {
 598         fCurrentChannel = iChannel;
 599         fCurrentSector  = icsector;
 600         fCurrentRow     = icRow;
 601     }
 602
 603     //process last pad if we change to a new one
 604     if ( iChannel != fCurrentChannel ){
 605         ProcessPad();
 606         fCurrentChannel = iChannel;
 607         fCurrentSector  = icsector;
 608         fCurrentRow     = icRow;
 609     }
 610
 611     //fill signals for current pad
 612     fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
 613     if ( csignal > fMaxPadSignal ){
 614         fMaxPadSignal = csignal;
 615         fMaxTimeBin   = icTimeBin;
 616     }
 617     return 0;
 618 }
 619 //_____________________________________________________________________
 620 void AliTPCCalibCE::FindPedestal(Float_t part)
 621 {
 622     //
 623     // find pedestal and noise for the current pad. Use either database or
 624     // truncated mean with part*100%
 625     //
 626     Bool_t noPedestal = kTRUE;
 627
 628     //use pedestal database if set
 629     if (fPedestalTPC&&fPadNoiseTPC){
 630         //only load new pedestals if the sector has changed
 631         if ( fCurrentSector!=fLastSector ){
 632             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 633             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 634             fLastSector=fCurrentSector;
 635         }
 636
 637         if ( fPedestalROC&&fPadNoiseROC ){
 638             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*fIsZeroSuppressed;
 639             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 640             noPedestal   = kFALSE;
 641         }
 642
 643     }
 644
 645     //if we are not running with pedestal database, or for the current sector there is no information
 646     //available, calculate the pedestal and noise on the fly
 647     if ( noPedestal ) {
 648         fPadPedestal = 0;
 649         fPadNoise    = 0;
 650         if ( fIsZeroSuppressed ) return;
 651         const Int_t kPedMax = 100;  //maximum pedestal value
 652         Float_t  max    =  0;
 653         Float_t  maxPos =  0;
 654         Int_t    median =  -1;
 655         Int_t    count0 =  0;
 656         Int_t    count1 =  0;
 657         //
 658         Float_t padSignal=0;
 659         //
 660         UShort_t histo[kPedMax];
 661         memset(histo,0,kPedMax*sizeof(UShort_t));
 662
 663         //fill pedestal histogram
 664         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 665             padSignal = fPadSignal.GetMatrixArray()[i];
 666             if (padSignal<=0) continue;
 667             if (padSignal>max && i>10) {
 668                 max = padSignal;
 669                 maxPos = i;
 670             }
 671             if (padSignal>kPedMax-1) continue;
 672             histo[int(padSignal+0.5)]++;
 673             count0++;
 674         }
 675         //find median
 676         for (Int_t i=1; i<kPedMax; ++i){
 677             if (count1<count0*0.5) median=i;
 678             count1+=histo[i];
 679         }
 680         // truncated mean
 681         //
 682         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 683         //
 684         for (Int_t idelta=1; idelta<10; ++idelta){
 685             if (median-idelta<=0) continue;
 686             if (median+idelta>kPedMax) continue;
 687             if (count<part*count1){
 688                 count+=histo[median-idelta];
 689                 mean +=histo[median-idelta]*(median-idelta);
 690                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 691                 count+=histo[median+idelta];
 692                 mean +=histo[median+idelta]*(median+idelta);
 693                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 694             }
 695         }
 696         if ( count > 0 ) {
 697             mean/=count;
 698             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 699             fPadPedestal = mean;
 700             fPadNoise    = rms;
 701         }
 702     }
 703 }
 704 //_____________________________________________________________________
 705 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
 706 {
 707     //
 708     //  Find position, signal width and height of the CE signal (last signal)
 709     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 710     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 711     //
 712
 713     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 714     Int_t   cemaxpos       = 0;
 715     Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise;  // threshold for the signal sum
 716     const Int_t    kCemin  = 4;             // range for the analysis of the ce signal +- channels from the peak
 717     const Int_t    kCemax  = 7;
 718
 719     Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
 720
 721     // find maximum closest to the sector mean from the last event
 722     for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
 723         // get sector mean of last event
 724         Float_t tmean = (*GetTMeanEvents(fCurrentSector))[fNevents-1];
 725             if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
 726                 minDist  = tmean-maxima[imax];
 727                 cemaxpos = (Int_t)maxima[imax];
 728             }
 729     }
 730
 731     if (cemaxpos!=0){
 732         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 733         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
 734             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
 735                 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 736                 if (signal>0) {
 737                     ceTime+=signal*(i+0.5);
 738                     ceRMS +=signal*(i+0.5)*(i+0.5);
 739                     ceQsum+=signal;
 740                 }
 741             }
 742         }
 743     }
 744     if (ceQmax&&ceQsum>ceSumThreshold) {
 745         ceTime/=ceQsum;
 746         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 747         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 748         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 749
 750         //Normalise Q to pad area of irocs
 751         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
 752
 753         ceQsum/=norm;
 754         fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
 755         fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
 756     } else {
 757         ceQmax=0;
 758         ceTime=0;
 759         ceRMS =0;
 760         ceQsum=0;
 761     }
 762     param[0] = ceQmax;
 763     param[1] = ceTime;
 764     param[2] = ceRMS;
 765     qSum     = ceQsum;
 766 }
 767 //_____________________________________________________________________
 768 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
 769 {
 770     //
 771     // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
 772     // and 'tplus' timebins after 'pos'
 773     //
 774     if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
 775     for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
 776         if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
 777     for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
 778         if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
 779             ++iTime2;
 780         if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
 781     }
 782     return kTRUE;
 783 }
 784 //_____________________________________________________________________
 785 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
 786 {
 787     //
 788     // Find local maxima on the pad signal and Histogram them
 789     //
 790   Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal
 791     Int_t   count       = 0;
 792 //    Int_t   tminus      = 2;
 793 //    Int_t   tplus       = 3;
 794     for (Int_t i=fLastTimeBin-fPeakPlus-1; i>=fFirstTimeBin+fPeakMinus; --i){
 795         if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,fPeakMinus,fPeakPlus) ){
 796           if (count<maxima.GetNrows()){
 797             maxima.GetMatrixArray()[count++]=i;
 798             GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
 799           }
 800         }
 801     }
 802 }
 803 //_____________________________________________________________________
 804 void AliTPCCalibCE::ProcessPad()
 805 {
 806     //
 807     //  Process data of current pad
 808     //
 809     FindPedestal();
 810
 811     TVectorF maxima(15);     // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
 812                              // + central electrode and possibly post peaks from the CE signal
 813                              // however if we are on a high noise pad a lot more peaks due to the noise might occur
 814     FindLocalMaxima(maxima);
 815     if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
 816
 817     if ( !GetTMeanEvents(fCurrentSector) ) return; //return if we don't have time 0 info, eg if only one side has laser
 818
 819     TVectorD param(3);
 820     Float_t  qSum;
 821     FindCESignal(param, qSum, maxima);
 822
 823     Double_t meanT  = param[1];
 824     Double_t sigmaT = param[2];
 825
 826     //Fill Event T0 counter
 827     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 828
 829     //Fill Q histogram
 830     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
 831
 832     //Fill RMS histogram
 833     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 834
 835
 836     //Fill debugging info
 837     if ( fDebugLevel>0 ){
 838         (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
 839         (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
 840         (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
 841     }
 842
 843     ResetPad();
 844 }
 845 //_____________________________________________________________________
 846 void AliTPCCalibCE::EndEvent()
 847 {
 848     //
 849     //  Process data of current pad
 850     //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
 851     //  before the EndEvent function to set the event timestamp and number!!!
 852     //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
 853     //  function was called
 854     //
 855
 856     //check if last pad has allready been processed, if not do so
 857     if ( fMaxTimeBin>-1 ) ProcessPad();
 858
 859 //    AliDebug(5,
 860
 861     TVectorD param(3);
 862     TMatrixD dummy(3,3);
 863 //    TVectorF vMeanTime(72);
 864 //    TVectorF vMeanQ(72);
 865     AliTPCCalROC *calIroc=new AliTPCCalROC(0);
 866     AliTPCCalROC *calOroc=new AliTPCCalROC(36);
 867
 868     //find mean time0 offset for side A and C
 869     Double_t time0Side[2];       //time0 for side A:0 and C:1
 870     Double_t time0SideCount[2];  //time0 counter for side A:0 and C:1
 871     time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
 872     for ( Int_t iSec = 0; iSec<72; ++iSec ){
 873         time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
 874         time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
 875     }
 876     if ( time0SideCount[0] >0  )
 877         time0Side[0]/=time0SideCount[0];
 878     if ( time0SideCount[1] >0 )
 879         time0Side[1]/=time0SideCount[1];
 880     // end find time0 offset
 881
 882     Int_t nSecMeanT=0;
 883     //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
 884     for ( Int_t iSec = 0; iSec<72; ++iSec ){
 885       //find median and then calculate the mean around it
 886         TH1S *hMeanT    = GetHistoTmean(iSec); //histogram with local maxima position information
 887         if ( !hMeanT ) continue;
 888         //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
 889         if ( hMeanT->GetEntries() < fROC->GetNChannels(iSec)*2/3 ){
 890             hMeanT->Reset();
 891             AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
 892           continue;
 893         }
 894
 895         Double_t entries = hMeanT->GetEntries();
 896         Double_t sum     = 0;
 897         Short_t *arr     = hMeanT->GetArray()+1;
 898         Int_t ibin=0;
 899         for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
 900             sum+=arr[ibin];
 901             if ( sum>=(entries/2.) ) break;
 902         }
 903         Int_t delta = 4;
 904         Int_t firstBin = fFirstTimeBin+ibin-delta;
 905         Int_t lastBin  = fFirstTimeBin+ibin+delta;
 906         if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
 907         if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
 908         Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
 909
 910         // check boundaries for ebye info of mean time
 911         TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
 912         Int_t vSize=vMeanTime->GetNrows();
 913         if ( vSize < fNevents+1 )
 914             vMeanTime->ResizeTo(vSize+100);
 915
 916         vMeanTime->GetMatrixArray()[fNevents]=median;
 917         nSecMeanT++;
 918       // end find median
 919
 920         TVectorF *vTimes = GetPadTimesEvent(iSec);
 921         if ( !vTimes ) continue;                     //continue if no time information for this sector is available
 922
 923
 924         AliTPCCalROC calIrocOutliers(0);
 925         AliTPCCalROC calOrocOutliers(36);
 926
 927         // calculate mean Q of the sector
 928         Float_t meanQ = 0;
 929         if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
 930         TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
 931         if ( vSize < fNevents+1 )           // vSize is the same as for vMeanTime!
 932             vMeanQ->ResizeTo(vSize+100);
 933
 934         vMeanQ->GetMatrixArray()[fNevents]=meanQ;
 935
 936         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
 937             Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
 938
 939             //set values for temporary roc calibration class
 940             if ( iSec < 36 ) {
 941                 calIroc->SetValue(iChannel, time);
 942                 if ( time == 0 ) calIrocOutliers.SetValue(iChannel,1);
 943
 944             } else {
 945                 calOroc->SetValue(iChannel, time);
 946                 if ( time == 0 ) calOrocOutliers.SetValue(iChannel,1);
 947             }
 948
 949             if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
 950                 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
 951
 952
 953
 954             //-------------------------------  Debug start  ------------------------------
 955             if ( fDebugLevel>0 ){
 956                 if ( !fDebugStreamer ) {
 957                         //debug stream
 958                     TDirectory *backup = gDirectory;
 959                     fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
 960                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 961                 }
 962
 963                 Int_t row=0;
 964                 Int_t pad=0;
 965                 Int_t padc=0;
 966
 967                 Float_t q   = (*GetPadQEvent(iSec))[iChannel];
 968                 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
 969
 970                 UInt_t channel=iChannel;
 971                 Int_t sector=iSec;
 972
 973                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 974                 pad = channel-fROC->GetRowIndexes(sector)[row];
 975                 padc = pad-(fROC->GetNPads(sector,row)/2);
 976
 977 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
 978 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
 979 //                                  fLastTimeBin-fFirstTimeBin,
 980 //                                  fFirstTimeBin,fLastTimeBin);
 981 //              h1->SetDirectory(0);
 982 //
 983 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 984 //                  h1->Fill(i,fPadSignal(i));
 985
 986                 Double_t t0Sec = 0;
 987                 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
 988                     t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
 989                 Double_t t0Side = time0Side[(iSec/18)%2];
 990                 (*fDebugStreamer) << "DataPad" <<
 991                     "Event=" << fNevents <<
 992                     "RunNumber=" << fRunNumber <<
 993                     "TimeStamp="   << fTimeStamp <<
 994                     "Sector="<< sector <<
 995                     "Row="   << row<<
 996                     "Pad="   << pad <<
 997                     "PadC="  << padc <<
 998                     "PadSec="<< channel <<
 999                     "Time0Sec="  << t0Sec <<
1000                     "Time0Side=" << t0Side <<
1001                     "Time="  << time <<
1002                     "RMS="   << rms <<
1003                     "Sum="   << q <<
1004                     "MeanQ=" << meanQ <<
1005                     //              "hist.=" << h1 <<
1006                     "\n";
1007
1008                 //              delete h1;
1009
1010             }
1011             //-----------------------------  Debug end  ------------------------------
1012         }// end channel loop
1013
1014         TVectorD paramPol1(3);
1015         TVectorD paramPol2(6);
1016         TMatrixD matPol1(3,3);
1017         TMatrixD matPol2(6,6);
1018         Float_t  chi2Pol1=0;
1019         Float_t  chi2Pol2=0;
1020
1021         if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1022             if ( iSec < 36 ){
1023                 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1024                 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1025             } else {
1026                 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1027                 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1028             }
1029
1030             GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1031             GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1032         }
1033
1034         //-------------------------------  Debug start  ------------------------------
1035         if ( fDebugLevel>0 ){
1036             if ( !fDebugStreamer ) {
1037                 //debug stream
1038                 TDirectory *backup = gDirectory;
1039                 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1040                 if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1041             }
1042             (*fDebugStreamer) << "DataRoc" <<
1043 //              "Event=" << fEvent <<
1044                 "RunNumber=" << fRunNumber <<
1045                 "TimeStamp="   << fTimeStamp <<
1046                 "Sector="<< iSec <<
1047                 "hMeanT.=" << hMeanT <<
1048                 "median=" << median <<
1049                 "paramPol1.=" << &paramPol1 <<
1050                 "paramPol2.=" << &paramPol2 <<
1051                 "matPol1.="   << &matPol1 <<
1052                 "matPol2.="   << &matPol2 <<
1053                 "chi2Pol1="   << chi2Pol1 <<
1054                 "chi2Pol2="   << chi2Pol2 <<
1055                 "\n";
1056         }
1057         //-------------------------------  Debug end  ------------------------------
1058         hMeanT->Reset();
1059     }// end sector loop
1060     //return if no sector has a valid mean time
1061     if ( nSecMeanT == 0 ) return;
1062
1063
1064 //    fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1065 //    fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1066     if ( fVEventTime.GetNrows() < fNevents+1 ) {
1067         fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1068         fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1069     }
1070     fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1071     fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1072
1073     fNevents++;
1074     fOldRunNumber = fRunNumber;
1075
1076     delete calIroc;
1077     delete calOroc;
1078 }
1079 //_____________________________________________________________________
1080 Bool_t AliTPCCalibCE::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
1081 {
1082   //
1083   // Event Processing loop - AliTPCRawStreamFast
1084   //
1085   ResetEvent();
1086   Bool_t withInput = kFALSE;
1087   while ( rawStreamFast->NextDDL() ){
1088     while ( rawStreamFast->NextChannel() ){
1089       Int_t isector  = rawStreamFast->GetSector();                       //  current sector
1090       Int_t iRow     = rawStreamFast->GetRow();                          //  current row
1091       Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
1092
1093       while ( rawStreamFast->NextBunch() ){
1094         Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
1095         Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
1096         for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
1097           Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
1098           Update(isector,iRow,iPad,iTimeBin+1,signal);
1099           withInput = kTRUE;
1100         }
1101       }
1102     }
1103   }
1104   if (withInput){
1105     EndEvent();
1106   }
1107   return withInput;
1108 }
1109 //_____________________________________________________________________
1110 Bool_t AliTPCCalibCE::ProcessEventFast(AliRawReader *rawReader)
1111 {
1112   //
1113   //  Event processing loop using the fast raw stream algorithm- AliRawReader
1114   //
1115
1116   //printf("ProcessEventFast - raw reader\n");
1117
1118   AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1119   if (eventHeader){
1120       fTimeStamp   = eventHeader->Get("Timestamp");
1121       fRunNumber = eventHeader->Get("RunNb");
1122   }
1123   fEventId = *rawReader->GetEventId();
1124
1125   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
1126   Bool_t res=ProcessEventFast(rawStreamFast);
1127   delete rawStreamFast;
1128   return res;
1129
1130 }
1131 //_____________________________________________________________________
1132 Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream)
1133 {
1134   //
1135   // Event Processing loop - AliTPCRawStream
1136   // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
1137   //
1138
1139   ResetEvent();
1140
1141   Bool_t withInput = kFALSE;
1142
1143   while (rawStream->Next()) {
1144
1145       Int_t isector  = rawStream->GetSector();                       //  current sector
1146       Int_t iRow     = rawStream->GetRow();                          //  current row
1147       Int_t iPad     = rawStream->GetPad();                          //  current pad
1148       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
1149       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
1150
1151       Update(isector,iRow,iPad,iTimeBin,signal);
1152       withInput = kTRUE;
1153   }
1154
1155   if (withInput){
1156       EndEvent();
1157   }
1158
1159   return withInput;
1160 }
1161 //_____________________________________________________________________
1162 Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
1163 {
1164   //
1165   //  Event processing loop - AliRawReader
1166   //
1167
1168
1169   AliTPCRawStream rawStream(rawReader,(AliAltroMapping**)fMapping);
1170     AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1171     if (eventHeader){
1172         fTimeStamp   = eventHeader->Get("Timestamp");
1173         fRunNumber = eventHeader->Get("RunNb");
1174     }
1175     fEventId = *rawReader->GetEventId();
1176
1177
1178     rawReader->Select("TPC");
1179
1180   return ProcessEvent(&rawStream);
1181 }
1182 //_____________________________________________________________________
1183 Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
1184 {
1185   //
1186   //  Event processing loop - date event
1187   //
1188     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1189     Bool_t result=ProcessEvent(rawReader);
1190     delete rawReader;
1191     return result;
1192
1193 }
1194 //_____________________________________________________________________
1195 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1196                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
1197                                   Char_t *type, Bool_t force)
1198 {
1199     //
1200     // return pointer to TH2S histogram of 'type'
1201     // if force is true create a new histogram if it doesn't exist allready
1202     //
1203     if ( !force || arr->UncheckedAt(sector) )
1204         return (TH2S*)arr->UncheckedAt(sector);
1205
1206     // if we are forced and histogram doesn't exist yet create it
1207     Char_t name[255], title[255];
1208
1209     sprintf(name,"hCalib%s%.2d",type,sector);
1210     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1211
1212     // new histogram with Q calib information. One value for each pad!
1213     TH2S* hist = new TH2S(name,title,
1214                           nbinsY, ymin, ymax,
1215                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1216     hist->SetDirectory(0);
1217     arr->AddAt(hist,sector);
1218     return hist;
1219 }
1220 //_____________________________________________________________________
1221 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1222 {
1223     //
1224     // return pointer to T0 histogram
1225     // if force is true create a new histogram if it doesn't exist allready
1226     //
1227     TObjArray *arr = &fHistoT0Array;
1228     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1229 }
1230 //_____________________________________________________________________
1231 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1232 {
1233     //
1234     // return pointer to Q histogram
1235     // if force is true create a new histogram if it doesn't exist allready
1236     //
1237     TObjArray *arr = &fHistoQArray;
1238     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1239 }
1240 //_____________________________________________________________________
1241 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1242 {
1243     //
1244     // return pointer to Q histogram
1245     // if force is true create a new histogram if it doesn't exist allready
1246     //
1247     TObjArray *arr = &fHistoRMSArray;
1248     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1249 }
1250 //_____________________________________________________________________
1251 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1252                               Char_t *type, Bool_t force)
1253 {
1254     //
1255     // return pointer to TH1S histogram
1256     // if force is true create a new histogram if it doesn't exist allready
1257     //
1258     if ( !force || arr->UncheckedAt(sector) )
1259         return (TH1S*)arr->UncheckedAt(sector);
1260
1261     // if we are forced and histogram doesn't yes exist create it
1262     Char_t name[255], title[255];
1263
1264     sprintf(name,"hCalib%s%.2d",type,sector);
1265     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1266
1267     // new histogram with calib information. One value for each pad!
1268     TH1S* hist = new TH1S(name,title,
1269                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1270     hist->SetDirectory(0);
1271     arr->AddAt(hist,sector);
1272     return hist;
1273 }
1274 //_____________________________________________________________________
1275 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1276 {
1277     //
1278     // return pointer to Q histogram
1279     // if force is true create a new histogram if it doesn't exist allready
1280     //
1281     TObjArray *arr = &fHistoTmean;
1282     return GetHisto(sector, arr, "LastTmean", force);
1283 }
1284 //_____________________________________________________________________
1285 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1286 {
1287     //
1288     // return pointer to Pad Info from 'arr' for the current event and sector
1289     // if force is true create it if it doesn't exist allready
1290     //
1291     if ( !force || arr->UncheckedAt(sector) )
1292         return (TVectorF*)arr->UncheckedAt(sector);
1293
1294     TVectorF *vect = new TVectorF(size);
1295     arr->AddAt(vect,sector);
1296     return vect;
1297 }
1298 //_____________________________________________________________________
1299 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1300 {
1301     //
1302     // return pointer to Pad Times Array for the current event and sector
1303     // if force is true create it if it doesn't exist allready
1304     //
1305     TObjArray *arr = &fPadTimesArrayEvent;
1306     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1307 }
1308 //_____________________________________________________________________
1309 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1310 {
1311     //
1312     // return pointer to Pad Q Array for the current event and sector
1313     // if force is true create it if it doesn't exist allready
1314     // for debugging purposes only
1315     //
1316
1317     TObjArray *arr = &fPadQArrayEvent;
1318     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1319 }
1320 //_____________________________________________________________________
1321 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1322 {
1323     //
1324     // return pointer to Pad RMS Array for the current event and sector
1325     // if force is true create it if it doesn't exist allready
1326     // for debugging purposes only
1327     //
1328     TObjArray *arr = &fPadRMSArrayEvent;
1329     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1330 }
1331 //_____________________________________________________________________
1332 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1333 {
1334     //
1335     // return pointer to Pad RMS Array for the current event and sector
1336     // if force is true create it if it doesn't exist allready
1337     // for debugging purposes only
1338     //
1339     TObjArray *arr = &fPadPedestalArrayEvent;
1340     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1341 }
1342 //_____________________________________________________________________
1343 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1344 {
1345     //
1346     // return pointer to the EbyE info of the mean arrival time for 'sector'
1347     // if force is true create it if it doesn't exist allready
1348     //
1349     TObjArray *arr = &fTMeanArrayEvent;
1350     return GetVectSector(sector,arr,100,force);
1351 }
1352 //_____________________________________________________________________
1353 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1354 {
1355     //
1356     // return pointer to the EbyE info of the mean arrival time for 'sector'
1357     // if force is true create it if it doesn't exist allready
1358     //
1359     TObjArray *arr = &fQMeanArrayEvent;
1360     return GetVectSector(sector,arr,100,force);
1361 }
1362 //_____________________________________________________________________
1363 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1364 {
1365     //
1366     // return pointer to ROC Calibration
1367     // if force is true create a new histogram if it doesn't exist allready
1368     //
1369     if ( !force || arr->UncheckedAt(sector) )
1370         return (AliTPCCalROC*)arr->UncheckedAt(sector);
1371
1372     // if we are forced and histogram doesn't yes exist create it
1373
1374     // new AliTPCCalROC for T0 information. One value for each pad!
1375     AliTPCCalROC *croc = new AliTPCCalROC(sector);
1376     arr->AddAt(croc,sector);
1377     return croc;
1378 }
1379 //_____________________________________________________________________
1380 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1381 {
1382     //
1383     // return pointer to Time 0 ROC Calibration
1384     // if force is true create a new histogram if it doesn't exist allready
1385     //
1386     TObjArray *arr = &fCalRocArrayT0;
1387     return GetCalRoc(sector, arr, force);
1388 }
1389 //_____________________________________________________________________
1390 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1391 {
1392     //
1393     // return pointer to the error of Time 0 ROC Calibration
1394     // if force is true create a new histogram if it doesn't exist allready
1395     //
1396     TObjArray *arr = &fCalRocArrayT0Err;
1397     return GetCalRoc(sector, arr, force);
1398 }
1399 //_____________________________________________________________________
1400 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1401 {
1402     //
1403     // return pointer to T0 ROC Calibration
1404     // if force is true create a new histogram if it doesn't exist allready
1405     //
1406     TObjArray *arr = &fCalRocArrayQ;
1407     return GetCalRoc(sector, arr, force);
1408 }
1409 //_____________________________________________________________________
1410 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1411 {
1412     //
1413     // return pointer to signal width ROC Calibration
1414     // if force is true create a new histogram if it doesn't exist allready
1415     //
1416     TObjArray *arr = &fCalRocArrayRMS;
1417     return GetCalRoc(sector, arr, force);
1418 }
1419 //_____________________________________________________________________
1420 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1421 {
1422     //
1423     // return pointer to Outliers
1424     // if force is true create a new histogram if it doesn't exist allready
1425     //
1426     TObjArray *arr = &fCalRocArrayOutliers;
1427     return GetCalRoc(sector, arr, force);
1428 }
1429 //_____________________________________________________________________
1430 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1431 {
1432     //
1433     // return pointer to TObjArray of fit parameters
1434     // if force is true create a new histogram if it doesn't exist allready
1435     //
1436     if ( !force || arr->UncheckedAt(sector) )
1437         return (TObjArray*)arr->UncheckedAt(sector);
1438
1439     // if we are forced and array doesn't yes exist create it
1440
1441     // new TObjArray for parameters
1442     TObjArray *newArr = new TObjArray;
1443     arr->AddAt(newArr,sector);
1444     return newArr;
1445 }
1446 //_____________________________________________________________________
1447 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1448 {
1449     //
1450     // return pointer to TObjArray of fit parameters from plane fit
1451     // if force is true create a new histogram if it doesn't exist allready
1452     //
1453     TObjArray *arr = &fParamArrayEventPol1;
1454     return GetParamArray(sector, arr, force);
1455 }
1456 //_____________________________________________________________________
1457 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1458 {
1459     //
1460     // return pointer to TObjArray of fit parameters from parabola fit
1461     // if force is true create a new histogram if it doesn't exist allready
1462     //
1463     TObjArray *arr = &fParamArrayEventPol2;
1464     return GetParamArray(sector, arr, force);
1465 }
1466 //_____________________________________________________________________
1467 void AliTPCCalibCE::ResetEvent()
1468 {
1469     //
1470     //  Reset global counters  -- Should be called before each event is processed
1471     //
1472     fLastSector=-1;
1473     fCurrentSector=-1;
1474     fCurrentRow=-1;
1475     fCurrentChannel=-1;
1476
1477     ResetPad();
1478
1479     fPadTimesArrayEvent.Delete();
1480     fPadQArrayEvent.Delete();
1481     fPadRMSArrayEvent.Delete();
1482     fPadPedestalArrayEvent.Delete();
1483
1484     for ( Int_t i=0; i<72; ++i ){
1485         fVTime0Offset.GetMatrixArray()[i]=0;
1486         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1487         fVMeanQ.GetMatrixArray()[i]=0;
1488         fVMeanQCounter.GetMatrixArray()[i]=0;
1489     }
1490 }
1491 //_____________________________________________________________________
1492 void AliTPCCalibCE::ResetPad()
1493 {
1494     //
1495     //  Reset pad infos -- Should be called after a pad has been processed
1496     //
1497     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1498         fPadSignal.GetMatrixArray()[i] = 0;
1499     fMaxTimeBin   = -1;
1500     fMaxPadSignal = -1;
1501     fPadPedestal  = -1;
1502     fPadNoise     = -1;
1503 }
1504 //_____________________________________________________________________
1505 void AliTPCCalibCE::Merge(AliTPCCalibCE *ce)
1506 {
1507     //
1508     //  Merge ce to the current AliTPCCalibCE
1509     //
1510
1511     //merge histograms
1512     for (Int_t iSec=0; iSec<72; ++iSec){
1513         TH2S *hRefQmerge   = ce->GetHistoQ(iSec);
1514         TH2S *hRefT0merge  = ce->GetHistoT0(iSec);
1515         TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1516
1517
1518         if ( hRefQmerge ){
1519             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1520             TH2S *hRefQ   = GetHistoQ(iSec);
1521             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1522             else {
1523                 TH2S *hist = new TH2S(*hRefQmerge);
1524                 hist->SetDirectory(0);
1525                 fHistoQArray.AddAt(hist, iSec);
1526             }
1527             hRefQmerge->SetDirectory(dir);
1528         }
1529         if ( hRefT0merge ){
1530             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1531             TH2S *hRefT0  = GetHistoT0(iSec);
1532             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1533             else {
1534                 TH2S *hist = new TH2S(*hRefT0merge);
1535                 hist->SetDirectory(0);
1536                 fHistoT0Array.AddAt(hist, iSec);
1537             }
1538             hRefT0merge->SetDirectory(dir);
1539         }
1540         if ( hRefRMSmerge ){
1541             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1542             TH2S *hRefRMS = GetHistoRMS(iSec);
1543             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1544             else {
1545                 TH2S *hist = new TH2S(*hRefRMSmerge);
1546                 hist->SetDirectory(0);
1547                 fHistoRMSArray.AddAt(hist, iSec);
1548             }
1549             hRefRMSmerge->SetDirectory(dir);
1550         }
1551
1552     }
1553
1554     // merge time information
1555
1556
1557     Int_t nCEevents = ce->GetNeventsProcessed();
1558     for (Int_t iSec=0; iSec<72; ++iSec){
1559         TObjArray *arrPol1CE  = ce->GetParamArrayPol1(iSec);
1560         TObjArray *arrPol2CE  = ce->GetParamArrayPol2(iSec);
1561         TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1562         TVectorF *vMeanQCE    = ce->GetQMeanEvents(iSec);
1563
1564         TObjArray *arrPol1  = 0x0;
1565         TObjArray *arrPol2  = 0x0;
1566         TVectorF *vMeanTime = 0x0;
1567         TVectorF *vMeanQ    = 0x0;
1568
1569         //resize arrays
1570         if ( arrPol1CE && arrPol2CE ){
1571             arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1572             arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1573             arrPol1->Expand(fNevents+nCEevents);
1574             arrPol2->Expand(fNevents+nCEevents);
1575         }
1576         if ( vMeanTimeCE && vMeanQCE ){
1577             vMeanTime = GetTMeanEvents(iSec,kTRUE);
1578             vMeanQ    = GetQMeanEvents(iSec,kTRUE);
1579             vMeanTime->ResizeTo(fNevents+nCEevents);
1580             vMeanQ->ResizeTo(fNevents+nCEevents);
1581         }
1582
1583
1584         for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1585             if ( arrPol1CE && arrPol2CE ){
1586                 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1587                 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1588                 if ( paramPol1 && paramPol2 ){
1589                     GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1590                     GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1591                 }
1592             }
1593             if ( vMeanTimeCE && vMeanQCE ){
1594                 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1595                 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1596             }
1597         }
1598     }
1599
1600
1601
1602     TVectorD*  eventTimes  = ce->GetEventTimes();
1603     TVectorD*  eventIds  = ce->GetEventIds();
1604     fVEventTime.ResizeTo(fNevents+nCEevents);
1605     fVEventNumber.ResizeTo(fNevents+nCEevents);
1606
1607     for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1608         Double_t evTime     = eventTimes->GetMatrixArray()[iEvent];
1609         Double_t evId       = eventIds->GetMatrixArray()[iEvent];
1610         fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1611         fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1612     }
1613     fNevents+=nCEevents; //increase event counter
1614
1615 }
1616 //_____________________________________________________________________
1617 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1618 {
1619     //
1620     // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1621     // xVariable:    0-event time, 1-event id, 2-internal event counter
1622     // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1623     // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1624     //                               0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1625     //                               not used for mean time and mean Q )
1626     // for an example see class description at the beginning
1627     //
1628
1629     Double_t *x = new Double_t[fNevents];
1630     Double_t *y = new Double_t[fNevents];
1631
1632     TVectorD *xVar = 0x0;
1633     TObjArray *aType = 0x0;
1634     Int_t npoints=0;
1635
1636     // sanity checks
1637     if ( (sector<0) || (sector>71) )      return 0x0;
1638     if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1639     if ( (fitType<0) || (fitType>3) )     return 0x0;
1640     if ( fitType==0 ){
1641         if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1642         aType = &fParamArrayEventPol1;
1643         if ( aType->At(sector)==0x0 ) return 0x0;
1644     }
1645     else if ( fitType==1 ){
1646         if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1647         aType = &fParamArrayEventPol2;
1648         if ( aType->At(sector)==0x0 ) return 0x0;
1649     }
1650
1651
1652     if ( xVariable == 0 ) xVar = &fVEventTime;
1653     if ( xVariable == 1 ) xVar = &fVEventNumber;
1654     if ( xVariable == 2 ) {
1655         xVar = new TVectorD(fNevents);
1656         for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1657     }
1658
1659     for (Int_t ievent =0; ievent<fNevents; ++ievent){
1660         if ( fitType<2 ){
1661             TObjArray *events = (TObjArray*)(aType->At(sector));
1662             if ( events->GetSize()<=ievent ) break;
1663             TVectorD *v = (TVectorD*)(events->At(ievent));
1664             if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1665         } else if (fitType == 2) {
1666             Double_t xValue=(*xVar)[ievent];
1667             Double_t yValue=(*GetTMeanEvents(sector))[ievent];
1668             if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1669         }else if (fitType == 3) {
1670             Double_t xValue=(*xVar)[ievent];
1671             Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1672             if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1673         }
1674     }
1675
1676     TGraph *gr = new TGraph(npoints);
1677     //sort xVariable increasing
1678     Int_t    *sortIndex = new Int_t[npoints];
1679     TMath::Sort(npoints,x,sortIndex);
1680     for (Int_t i=0;i<npoints;++i){
1681         gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1682     }
1683
1684
1685     if ( xVariable == 2 ) delete xVar;
1686     delete x;
1687     delete y;
1688     delete sortIndex;
1689     return gr;
1690 }
1691 //_____________________________________________________________________
1692 void AliTPCCalibCE::Analyse()
1693 {
1694     //
1695     //  Calculate calibration constants
1696     //
1697
1698     TVectorD paramQ(3);
1699     TVectorD paramT0(3);
1700     TVectorD paramRMS(3);
1701     TMatrixD dummy(3,3);
1702
1703     Float_t channelCounter=0;
1704     fMeanT0rms=0;
1705     fMeanQrms=0;
1706     fMeanRMSrms=0;
1707
1708     for (Int_t iSec=0; iSec<72; ++iSec){
1709         TH2S *hT0 = GetHistoT0(iSec);
1710         if (!hT0 ) continue;
1711
1712         AliTPCCalROC *rocQ     = GetCalRocQ  (iSec,kTRUE);
1713         AliTPCCalROC *rocT0    = GetCalRocT0 (iSec,kTRUE);
1714         AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1715         AliTPCCalROC *rocRMS   = GetCalRocRMS(iSec,kTRUE);
1716         AliTPCCalROC *rocOut   = GetCalRocOutliers(iSec,kTRUE);
1717
1718         TH2S *hQ   = GetHistoQ(iSec);
1719         TH2S *hRMS = GetHistoRMS(iSec);
1720
1721         Short_t *arrayhQ   = hQ->GetArray();
1722         Short_t *arrayhT0  = hT0->GetArray();
1723         Short_t *arrayhRMS = hRMS->GetArray();
1724
1725         UInt_t nChannels = fROC->GetNChannels(iSec);
1726
1727         //debug
1728         Int_t row=0;
1729         Int_t pad=0;
1730         Int_t padc=0;
1731         //! debug
1732
1733         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1734
1735
1736             Float_t cogTime0 = -1000;
1737             Float_t cogQ     = -1000;
1738             Float_t cogRMS   = -1000;
1739             Float_t cogOut   = 0;
1740             Float_t rms      = 0;
1741             Float_t rmsT0    = 0;
1742
1743
1744             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1745             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1746             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1747
1748             cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
1749             fMeanQrms+=rms;
1750             cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
1751             fMeanT0rms+=rmsT0;
1752             cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
1753             fMeanRMSrms+=rms;
1754             channelCounter++;
1755
1756             /*
1757              //outlier specifications
1758             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1759                 cogOut = 1;
1760                 cogTime0 = 0;
1761                 cogQ     = 0;
1762                 cogRMS   = 0;
1763             }
1764 */
1765             rocQ->SetValue(iChannel, cogQ*cogQ);
1766             rocT0->SetValue(iChannel, cogTime0);
1767             rocT0Err->SetValue(iChannel, rmsT0);
1768             rocRMS->SetValue(iChannel, cogRMS);
1769             rocOut->SetValue(iChannel, cogOut);
1770
1771
1772             //debug
1773             if ( fDebugLevel > 0 ){
1774                 if ( !fDebugStreamer ) {
1775                         //debug stream
1776                     TDirectory *backup = gDirectory;
1777                     fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1778                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1779                 }
1780
1781                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1782                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1783                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1784
1785                 (*fDebugStreamer) << "DataEnd" <<
1786                     "Sector="  << iSec      <<
1787                     "Pad="     << pad       <<
1788                     "PadC="    << padc      <<
1789                     "Row="     << row       <<
1790                     "PadSec="  << iChannel   <<
1791                     "Q="       << cogQ      <<
1792                     "T0="      << cogTime0  <<
1793                     "RMS="     << cogRMS    <<
1794                     "\n";
1795             }
1796             //! debug
1797
1798         }
1799
1800     }
1801     if ( channelCounter>0 ){
1802         fMeanT0rms/=channelCounter;
1803         fMeanQrms/=channelCounter;
1804         fMeanRMSrms/=channelCounter;
1805     }
1806     if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
1807 //    delete fDebugStreamer;
1808 //    fDebugStreamer = 0x0;
1809 }
1810 //_____________________________________________________________________
1811 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1812 {
1813     //
1814     //  Write class to file
1815     //
1816
1817     TString sDir(dir);
1818     TString option;
1819
1820     if ( append )
1821         option = "update";
1822     else
1823         option = "recreate";
1824
1825     TDirectory *backup = gDirectory;
1826     TFile f(filename,option.Data());
1827     f.cd();
1828     if ( !sDir.IsNull() ){
1829         f.mkdir(sDir.Data());
1830         f.cd(sDir);
1831     }
1832     this->Write();
1833     f.Close();
1834
1835     if ( backup ) backup->cd();
1836 }