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