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