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