TPC/AliTPCCalibCE.cxx

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