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