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 <TH1.h>
 264 #include <TH1F.h>
 265 #include <TH2S.h>
 266 #include <TF1.h>
 267 #include <TString.h>
 268 #include <TVectorF.h>
 269 #include <TVectorD.h>
 270 #include <TVector3.h>
 271 #include <TMatrixD.h>
 272 #include <TMath.h>
 273 #include <TGraph.h>
 274 #include <TGraphErrors.h>
 275 #include <TString.h>
 276 #include <TMap.h>
 277 #include <TDirectory.h>
 278 #include <TSystem.h>
 279 #include <TFile.h>
 280 #include <TCollection.h>
 281 #include <TTimeStamp.h>
 282 #include <TList.h>
 283 #include <TKey.h>
 284
 285 //AliRoot includes
 286 #include "AliLog.h"
 287 #include "AliRawReader.h"
 288 #include "AliRawReaderRoot.h"
 289 #include "AliRawReaderDate.h"
 290 #include "AliRawEventHeaderBase.h"
 291 #include "AliTPCRawStream.h"
 292 #include "AliTPCRawStreamFast.h"
 293 #include "AliTPCCalROC.h"
 294 #include "AliTPCCalPad.h"
 295 #include "AliTPCROC.h"
 296 #include "AliTPCParam.h"
 297 #include "AliTPCCalibCE.h"
 298 #include "AliMathBase.h"
 299 #include "AliTPCTransform.h"
 300 #include "AliTPCLaserTrack.h"
 301 #include "TTreeStream.h"
 302
 303 #include "AliCDBManager.h"
 304 #include "AliCDBEntry.h"
 305 //date
 306 #include "event.h"
 307 ClassImp(AliTPCCalibCE)
 308
 309
 310 AliTPCCalibCE::AliTPCCalibCE() :
 311   AliTPCCalibRawBase(),
 312   fNbinsT0(200),
 313   fXminT0(-5),
 314   fXmaxT0(5),
 315   fNbinsQ(200),
 316   fXminQ(1),
 317   fXmaxQ(40),
 318   fNbinsRMS(100),
 319   fXminRMS(0.1),
 320   fXmaxRMS(5.1),
 321   fPeakDetMinus(2),
 322   fPeakDetPlus(3),
 323   fPeakIntMinus(2),
 324   fPeakIntPlus(2),
 325   fNoiseThresholdMax(5.),
 326   fNoiseThresholdSum(8.),
 327   fIsZeroSuppressed(kFALSE),
 328   fLastSector(-1),
 329   fSecRejectRatio(.4),
 330   fParam(new AliTPCParam),
 331   fPedestalTPC(0x0),
 332   fPadNoiseTPC(0x0),
 333   fPedestalROC(0x0),
 334   fPadNoiseROC(0x0),
 335   fCalRocArrayT0(72),
 336   fCalRocArrayT0Err(72),
 337   fCalRocArrayQ(72),
 338   fCalRocArrayRMS(72),
 339   fCalRocArrayOutliers(72),
 340   fHistoQArray(72),
 341   fHistoT0Array(72),
 342   fHistoRMSArray(72),
 343   fMeanT0rms(0),
 344   fMeanQrms(0),
 345   fMeanRMSrms(0),
 346   fHistoTmean(72),
 347   fParamArrayEventPol1(72),
 348   fParamArrayEventPol2(72),
 349   fTMeanArrayEvent(72),
 350   fQMeanArrayEvent(72),
 351   fVEventTime(1000),
 352   fVEventNumber(1000),
 353   fVTime0SideA(1000),
 354   fVTime0SideC(1000),
 355   fEventId(-1),
 356   fOldRunNumber(0),
 357   fPadTimesArrayEvent(72),
 358   fPadQArrayEvent(72),
 359   fPadRMSArrayEvent(72),
 360   fPadPedestalArrayEvent(72),
 361   fCurrentChannel(-1),
 362   fCurrentSector(-1),
 363   fCurrentRow(-1),
 364   fMaxPadSignal(-1),
 365   fMaxTimeBin(-1),
 366 //   fPadSignal(1024),
 367   fPadPedestal(0),
 368   fPadNoise(0),
 369   fVTime0Offset(72),
 370   fVTime0OffsetCounter(72),
 371   fVMeanQ(72),
 372   fVMeanQCounter(72),
 373   fCurrentCETimeRef(0),
 374   fProcessOld(kTRUE),
 375   fProcessNew(kFALSE),
 376   fAnalyseNew(kTRUE),
 377   fHnDrift(0x0),
 378   fArrHnDrift(100),
 379   fTimeBursts(100),
 380   fArrFitGraphs(0x0)
 381 {
 382   //
 383   // AliTPCSignal default constructor
 384   //
 385   SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
 386   fFirstTimeBin=650;
 387   fLastTimeBin=1030;
 388   fParam->Update();
 389   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 390   for (Int_t i=0;i<5;++i){
 391     fPeaks[i]=0;
 392     fPeakWidths[i]=0;
 393   }
 394   for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
 395 }
 396 //_____________________________________________________________________
 397 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
 398   AliTPCCalibRawBase(sig),
 399   fNbinsT0(sig.fNbinsT0),
 400   fXminT0(sig.fXminT0),
 401   fXmaxT0(sig.fXmaxT0),
 402   fNbinsQ(sig.fNbinsQ),
 403   fXminQ(sig.fXminQ),
 404   fXmaxQ(sig.fXmaxQ),
 405   fNbinsRMS(sig.fNbinsRMS),
 406   fXminRMS(sig.fXminRMS),
 407   fXmaxRMS(sig.fXmaxRMS),
 408   fPeakDetMinus(sig.fPeakDetMinus),
 409   fPeakDetPlus(sig.fPeakDetPlus),
 410   fPeakIntMinus(sig.fPeakIntMinus),
 411   fPeakIntPlus(sig.fPeakIntPlus),
 412   fNoiseThresholdMax(sig.fNoiseThresholdMax),
 413   fNoiseThresholdSum(sig.fNoiseThresholdSum),
 414   fIsZeroSuppressed(sig.fIsZeroSuppressed),
 415   fLastSector(-1),
 416   fSecRejectRatio(.4),
 417   fParam(new AliTPCParam),
 418   fPedestalTPC(0x0),
 419   fPadNoiseTPC(0x0),
 420   fPedestalROC(0x0),
 421   fPadNoiseROC(0x0),
 422   fCalRocArrayT0(72),
 423   fCalRocArrayT0Err(72),
 424   fCalRocArrayQ(72),
 425   fCalRocArrayRMS(72),
 426   fCalRocArrayOutliers(72),
 427   fHistoQArray(72),
 428   fHistoT0Array(72),
 429   fHistoRMSArray(72),
 430   fMeanT0rms(sig.fMeanT0rms),
 431   fMeanQrms(sig.fMeanQrms),
 432   fMeanRMSrms(sig.fMeanRMSrms),
 433   fHistoTmean(72),
 434   fParamArrayEventPol1(72),
 435   fParamArrayEventPol2(72),
 436   fTMeanArrayEvent(72),
 437   fQMeanArrayEvent(72),
 438   fVEventTime(sig.fVEventTime),
 439   fVEventNumber(sig.fVEventNumber),
 440   fVTime0SideA(sig.fVTime0SideA),
 441   fVTime0SideC(sig.fVTime0SideC),
 442   fEventId(-1),
 443   fOldRunNumber(0),
 444   fPadTimesArrayEvent(72),
 445   fPadQArrayEvent(72),
 446   fPadRMSArrayEvent(72),
 447   fPadPedestalArrayEvent(72),
 448   fCurrentChannel(-1),
 449   fCurrentSector(-1),
 450   fCurrentRow(-1),
 451   fMaxPadSignal(-1),
 452   fMaxTimeBin(-1),
 453 //   fPadSignal(1024),
 454   fPadPedestal(0),
 455   fPadNoise(0),
 456   fVTime0Offset(72),
 457   fVTime0OffsetCounter(72),
 458   fVMeanQ(72),
 459   fVMeanQCounter(72),
 460   fCurrentCETimeRef(0),
 461   fProcessOld(sig.fProcessOld),
 462   fProcessNew(sig.fProcessNew),
 463   fAnalyseNew(sig.fAnalyseNew),
 464   fHnDrift(0x0),
 465   fArrHnDrift(100),
 466   fTimeBursts(100),
 467   fArrFitGraphs(0x0)
 468 {
 469   //
 470   // AliTPCSignal copy constructor
 471   //
 472   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 473
 474   for (Int_t iSec = 0; iSec < 72; ++iSec){
 475     const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 476     const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 477     const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 478     const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 479
 480     const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 481     const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 482     const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 483
 484     if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 485     if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 486     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 487     if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 488
 489     if ( hQ != 0x0 ){
 490       TH2S *hNew = new TH2S(*hQ);
 491       hNew->SetDirectory(0);
 492       fHistoQArray.AddAt(hNew,iSec);
 493     }
 494     if ( hT0 != 0x0 ){
 495       TH2S *hNew = new TH2S(*hT0);
 496       hNew->SetDirectory(0);
 497       fHistoT0Array.AddAt(hNew,iSec);
 498     }
 499     if ( hRMS != 0x0 ){
 500       TH2S *hNew = new TH2S(*hRMS);
 501       hNew->SetDirectory(0);
 502       fHistoRMSArray.AddAt(hNew,iSec);
 503     }
 504   }
 505
 506   //copy fit parameters event by event
 507   TObjArray *arr=0x0;
 508   for (Int_t iSec=0; iSec<72; ++iSec){
 509     arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
 510     if ( arr ){
 511       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 512       fParamArrayEventPol1.AddAt(arrEvents, iSec);
 513       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 514         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 515           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 516     }
 517
 518     arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
 519     if ( arr ){
 520       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 521       fParamArrayEventPol2.AddAt(arrEvents, iSec);
 522       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 523         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 524           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 525     }
 526
 527     TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
 528     TVectorF *vMeanQ    = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
 529     if ( vMeanTime )
 530       fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
 531     if ( vMeanQ )
 532       fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
 533   }
 534
 535
 536   fVEventTime.ResizeTo(sig.fVEventTime);
 537   fVEventNumber.ResizeTo(sig.fVEventNumber);
 538   fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
 539   fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
 540
 541   fParam->Update();
 542
 543   for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
 544     TObject *o=sig.fArrHnDrift.UncheckedAt(i);
 545     if (o){
 546       TObject *newo=o->Clone("fHnDrift");
 547       fArrHnDrift.AddAt(newo,i);
 548       if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
 549     }
 550   }
 551
 552   for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
 553     fTimeBursts[i]=sig.fTimeBursts[i];
 554   }
 555
 556   for (Int_t i=0;i<5;++i){
 557     fPeaks[i]=sig.fPeaks[i];
 558     fPeakWidths[i]=sig.fPeakWidths[i];
 559   }
 560   if (sig.fArrFitGraphs) {
 561     fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
 562     fArrFitGraphs->SetOwner();
 563   }
 564
 565   for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
 566
 567 }
 568 //_____________________________________________________________________
 569 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
 570   AliTPCCalibRawBase(),
 571   fNbinsT0(200),
 572   fXminT0(-5),
 573   fXmaxT0(5),
 574   fNbinsQ(200),
 575   fXminQ(1),
 576   fXmaxQ(40),
 577   fNbinsRMS(100),
 578   fXminRMS(0.1),
 579   fXmaxRMS(5.1),
 580   fPeakDetMinus(2),
 581   fPeakDetPlus(3),
 582   fPeakIntMinus(2),
 583   fPeakIntPlus(2),
 584   fNoiseThresholdMax(5.),
 585   fNoiseThresholdSum(8.),
 586   fIsZeroSuppressed(kFALSE),
 587   fLastSector(-1),
 588   fSecRejectRatio(.4),
 589   fParam(new  AliTPCParam),
 590   fPedestalTPC(0x0),
 591   fPadNoiseTPC(0x0),
 592   fPedestalROC(0x0),
 593   fPadNoiseROC(0x0),
 594   fCalRocArrayT0(72),
 595   fCalRocArrayT0Err(72),
 596   fCalRocArrayQ(72),
 597   fCalRocArrayRMS(72),
 598   fCalRocArrayOutliers(72),
 599   fHistoQArray(72),
 600   fHistoT0Array(72),
 601   fHistoRMSArray(72),
 602   fMeanT0rms(0),
 603   fMeanQrms(0),
 604   fMeanRMSrms(0),
 605   fHistoTmean(72),
 606   fParamArrayEventPol1(72),
 607   fParamArrayEventPol2(72),
 608   fTMeanArrayEvent(72),
 609   fQMeanArrayEvent(72),
 610   fVEventTime(1000),
 611   fVEventNumber(1000),
 612   fVTime0SideA(1000),
 613   fVTime0SideC(1000),
 614   fEventId(-1),
 615   fOldRunNumber(0),
 616   fPadTimesArrayEvent(72),
 617   fPadQArrayEvent(72),
 618   fPadRMSArrayEvent(72),
 619   fPadPedestalArrayEvent(72),
 620   fCurrentChannel(-1),
 621   fCurrentSector(-1),
 622   fCurrentRow(-1),
 623   fMaxPadSignal(-1),
 624   fMaxTimeBin(-1),
 625 //   fPadSignal(1024),
 626   fPadPedestal(0),
 627   fPadNoise(0),
 628   fVTime0Offset(72),
 629   fVTime0OffsetCounter(72),
 630   fVMeanQ(72),
 631   fVMeanQCounter(72),
 632   fCurrentCETimeRef(0),
 633   fProcessOld(kTRUE),
 634   fProcessNew(kFALSE),
 635   fAnalyseNew(kTRUE),
 636   fHnDrift(0x0),
 637   fArrHnDrift(100),
 638   fTimeBursts(100),
 639   fArrFitGraphs(0x0)
 640 {
 641   //
 642   // constructor which uses a tmap as input to set some specific parameters
 643   //
 644   SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
 645   fFirstTimeBin=650;
 646   fLastTimeBin=1030;
 647   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 648   if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 649   if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
 650   if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
 651   if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
 652   if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
 653   if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
 654   if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
 655   if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
 656   if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
 657   if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
 658   if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
 659   if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
 660   if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
 661   if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
 662   if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
 663   if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
 664   if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
 665   if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
 666   if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
 667
 668   if (config->GetValue("ProcessOld")) fProcessOld = (Bool_t)((TObjString*)config->GetValue("ProcessOld"))->GetString().Atoi();
 669   if (config->GetValue("ProcessNew")) fProcessNew = (Bool_t)((TObjString*)config->GetValue("ProcessNew"))->GetString().Atoi();
 670   if (config->GetValue("AnalyseNew")) fAnalyseNew = (Bool_t)((TObjString*)config->GetValue("AnalyseNew"))->GetString().Atoi();
 671
 672   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 673   for (Int_t i=0;i<5;++i){
 674     fPeaks[i]=0;
 675     fPeakWidths[i]=0;
 676   }
 677
 678   fParam->Update();
 679   for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
 680 }
 681
 682 //_____________________________________________________________________
 683 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
 684 {
 685   //
 686   // assignment operator
 687   //
 688   if (&source == this) return *this;
 689   new (this) AliTPCCalibCE(source);
 690
 691   return *this;
 692 }
 693 //_____________________________________________________________________
 694 AliTPCCalibCE::~AliTPCCalibCE()
 695 {
 696   //
 697   // destructor
 698   //
 699
 700   fCalRocArrayT0.Delete();
 701   fCalRocArrayT0Err.Delete();
 702   fCalRocArrayQ.Delete();
 703   fCalRocArrayRMS.Delete();
 704   fCalRocArrayOutliers.Delete();
 705
 706   fHistoQArray.Delete();
 707   fHistoT0Array.Delete();
 708   fHistoRMSArray.Delete();
 709
 710   fHistoTmean.Delete();
 711
 712   fParamArrayEventPol1.Delete();
 713   fParamArrayEventPol2.Delete();
 714   fTMeanArrayEvent.Delete();
 715   fQMeanArrayEvent.Delete();
 716
 717   fPadTimesArrayEvent.Delete();
 718   fPadQArrayEvent.Delete();
 719   fPadRMSArrayEvent.Delete();
 720   fPadPedestalArrayEvent.Delete();
 721
 722   fArrHnDrift.SetOwner();
 723   fArrHnDrift.Delete();
 724
 725   if (fArrFitGraphs){
 726     fArrFitGraphs->SetOwner();
 727     delete fArrFitGraphs;
 728   }
 729 }
 730 //_____________________________________________________________________
 731 Int_t AliTPCCalibCE::Update(const Int_t icsector,
 732                                 const Int_t icRow,
 733                                 const Int_t icPad,
 734                                 const Int_t icTimeBin,
 735                                 const Float_t csignal)
 736 {
 737   //
 738   // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 739   // no extra analysis necessary. Assumes knowledge of the signal shape!
 740   // assumes that it is looped over consecutive time bins of one pad
 741   //
 742
 743   if (!fProcessOld) return 0;
 744   //temp
 745
 746   if (icRow<0) return 0;
 747   if (icPad<0) return 0;
 748   if (icTimeBin<0) return 0;
 749   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 750
 751   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 752
 753   //init first pad and sector in this event
 754   if ( fCurrentChannel == -1 ) {
 755     fLastSector=-1;
 756     fCurrentChannel = iChannel;
 757     fCurrentSector  = icsector;
 758     fCurrentRow     = icRow;
 759   }
 760
 761   //process last pad if we change to a new one
 762   if ( iChannel != fCurrentChannel ){
 763     ProcessPad();
 764     fLastSector=fCurrentSector;
 765     fCurrentChannel = iChannel;
 766     fCurrentSector  = icsector;
 767     fCurrentRow     = icRow;
 768   }
 769
 770   //fill signals for current pad
 771   fPadSignal[icTimeBin]=csignal;
 772   if ( csignal > fMaxPadSignal ){
 773     fMaxPadSignal = csignal;
 774     fMaxTimeBin   = icTimeBin;
 775   }
 776   return 0;
 777 }
 778
 779 //_____________________________________________________________________
 780 void AliTPCCalibCE::ProcessBunch(const Int_t sector, const Int_t row, const Int_t pad,
 781                   const Int_t length, const UInt_t startTimeBin, const UShort_t* signal)
 782 {
 783   //
 784   // new filling method to fill the THnSparse histogram
 785   //
 786
 787   //only in new processing mode
 788   if (!fProcessNew) return;
 789   //don't use the IROCs
 790   if (sector<36) return;
 791   //only bunches with reasonable length
 792   if (length<3||length>10) return;
 793
 794   UShort_t  timeBin    = (UShort_t)startTimeBin;
 795   //skip first laser layer
 796   if (timeBin<200) return;
 797   Double_t timeBurst=SetBurstHnDrift();
 798
 799   //after 1 event setup peak ranges
 800   if (fNevents==1&&fPeaks[4]==0) {
 801     // set time range
 802     fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
 803     FindLaserLayers();
 804     // set time range
 805     fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
 806   }
 807
 808   // After the first event only fill every 5th  bin in a row with the CE information
 809   Int_t padFill=pad;
 810   if (fPeaks[4]>100&&TMath::Abs((Short_t)fPeaks[4]-(Short_t)timeBin)<(Short_t)fPeakWidths[4]){
 811     Int_t mod=5;
 812     Int_t n=pad/mod;
 813     padFill=mod*n+mod/2;
 814   }
 815
 816   //noise removal
 817   if (!IsPeakInRange(timeBin+length/2)) return;
 818
 819   Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
 820       (Double_t)padFill,(Double_t)timeBin,timeBurst};
 821
 822   for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
 823     Float_t sig=(Float_t)signal[iTimeBin];
 824     x[3]=timeBin;
 825     fHnDrift->Fill(x,sig);
 826     --timeBin;
 827   }
 828 }
 829 //_____________________________________________________________________
 830 void AliTPCCalibCE::FindLaserLayers()
 831 {
 832   //
 833   // Find the laser layer positoins
 834   //
 835
 836
 837   //find CE signal position and width
 838   TH1D *hproj=fHnDrift->Projection(3);
 839   hproj->GetXaxis()->SetRangeUser(700,1030);
 840   Int_t maxbin=hproj->GetMaximumBin();
 841   Double_t binc=hproj->GetBinCenter(maxbin);
 842   hproj->GetXaxis()->SetRangeUser(binc-10,binc+10);
 843
 844   fPeaks[4]=(UShort_t)TMath::Nint(hproj->GetMean());
 845 //   fPeakWidths[4]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
 846   fPeakWidths[4]=(UShort_t)TMath::Nint(10.);
 847
 848   //other peaks
 849   Int_t timepos=fPeaks[4]-2*fPeakWidths[4];
 850   Int_t width=100;
 851
 852   for (Int_t i=3; i>=0; --i){
 853     hproj->GetXaxis()->SetRangeUser(timepos-width,timepos);
 854     fPeaks[i]=(UShort_t)TMath::Nint(hproj->GetMean());
 855     fPeakWidths[i]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
 856     width=250;
 857     timepos=fPeaks[i]-width/2;
 858   }
 859
 860   //check width and reset peak if >100
 861   for (Int_t i=0; i<5; ++i){
 862     if (fPeakWidths[i]>100) {
 863       fPeaks[i]=0;
 864       fPeakWidths[i]=0;
 865     }
 866   }
 867
 868   delete hproj;
 869 }
 870
 871 //_____________________________________________________________________
 872 void AliTPCCalibCE::FindPedestal(Float_t part)
 873 {
 874   //
 875     // find pedestal and noise for the current pad. Use either database or
 876     // truncated mean with part*100%
 877   //
 878   Bool_t noPedestal = kTRUE;
 879
 880     //use pedestal database if set
 881   if (fPedestalTPC&&fPadNoiseTPC){
 882         //only load new pedestals if the sector has changed
 883     if ( fCurrentSector!=fLastSector ){
 884       fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 885       fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 886     }
 887
 888     if ( fPedestalROC&&fPadNoiseROC ){
 889       fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
 890       fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 891       noPedestal   = kFALSE;
 892     }
 893
 894   }
 895
 896     //if we are not running with pedestal database, or for the current sector there is no information
 897     //available, calculate the pedestal and noise on the fly
 898   if ( noPedestal ) {
 899     fPadPedestal = 0;
 900     fPadNoise    = 0;
 901     if ( fIsZeroSuppressed ) return;
 902     const Int_t kPedMax = 100;  //maximum pedestal value
 903     Float_t  max    =  0;
 904     Float_t  maxPos =  0;
 905     Int_t    median =  -1;
 906     Int_t    count0 =  0;
 907     Int_t    count1 =  0;
 908     //
 909     Float_t padSignal=0;
 910     //
 911     UShort_t histo[kPedMax];
 912     memset(histo,0,kPedMax*sizeof(UShort_t));
 913
 914         //fill pedestal histogram
 915     for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 916       padSignal = fPadSignal[i];
 917       if (padSignal<=0) continue;
 918       if (padSignal>max && i>10) {
 919         max = padSignal;
 920         maxPos = i;
 921       }
 922       if (padSignal>kPedMax-1) continue;
 923       histo[int(padSignal+0.5)]++;
 924       count0++;
 925     }
 926         //find median
 927     for (Int_t i=1; i<kPedMax; ++i){
 928       if (count1<count0*0.5) median=i;
 929       count1+=histo[i];
 930     }
 931         // truncated mean
 932     //
 933     Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 934     //
 935     for (Int_t idelta=1; idelta<10; ++idelta){
 936       if (median-idelta<=0) continue;
 937       if (median+idelta>kPedMax) continue;
 938       if (count<part*count1){
 939         count+=histo[median-idelta];
 940         mean +=histo[median-idelta]*(median-idelta);
 941         rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 942         count+=histo[median+idelta];
 943         mean +=histo[median+idelta]*(median+idelta);
 944         rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 945       }
 946     }
 947     if ( count > 0 ) {
 948       mean/=count;
 949       rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 950       fPadPedestal = mean;
 951       fPadNoise    = rms;
 952     }
 953   }
 954 }
 955 //_____________________________________________________________________
 956 void AliTPCCalibCE::UpdateCETimeRef()
 957 {
 958   // Find the time reference of the last valid CE signal in sector
 959   // for irocs of the A-Side the reference of the corresponging OROC is returned
 960   // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
 961   if ( fLastSector == fCurrentSector ) return;
 962   Int_t sector=fCurrentSector;
 963   if ( sector < 18 ) sector+=36;
 964   fCurrentCETimeRef=0;
 965   TVectorF *vtRef = GetTMeanEvents(sector);
 966   if ( !vtRef ) return;
 967   Int_t vtRefSize= vtRef->GetNrows();
 968   if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
 969   else vtRefSize=fNevents;
 970   while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
 971   fCurrentCETimeRef=(*vtRef)[vtRefSize];
 972   AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
 973 }
 974 //_____________________________________________________________________
 975 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
 976 {
 977   //
 978     //  Find position, signal width and height of the CE signal (last signal)
 979     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 980     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 981   //
 982
 983   Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 984   Int_t   cemaxpos       = 0;
 985   Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise;  // threshold for the signal sum
 986   const Int_t    kCemin  = fPeakIntMinus;             // range for the analysis of the ce signal +- channels from the peak
 987   const Int_t    kCemax  = fPeakIntPlus;
 988
 989   Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
 990
 991     // find maximum closest to the sector mean from the last event
 992   for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
 993         // get sector mean of last event
 994     Float_t tmean = fCurrentCETimeRef;
 995     if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
 996       minDist  = tmean-maxima[imax];
 997       cemaxpos = (Int_t)maxima[imax];
 998     }
 999   }
1000 //   printf("L1 phase TB: %f\n",GetL1PhaseTB());
1001   if (cemaxpos!=0){
1002     ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
1003     for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
1004       if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
1005         Float_t signal = fPadSignal[i]-fPadPedestal;
1006         if (signal>0) {
1007           ceTime+=signal*(i+0.5);
1008           ceRMS +=signal*(i+0.5)*(i+0.5);
1009           ceQsum+=signal;
1010         }
1011       }
1012     }
1013   }
1014   if (ceQmax&&ceQsum>ceSumThreshold) {
1015     ceTime/=ceQsum;
1016     ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
1017     ceTime-=GetL1PhaseTB();
1018     fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
1019     fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
1020
1021   //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
1022   //                                the pick-up signal should scale with the pad area. In addition
1023   //                                the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
1024   //                                ratio 2/3. The pad area we express in cm2. We normalise the signal
1025   //                                to the OROC signal (factor 2/3 for the IROCs).
1026     Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
1027     if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
1028
1029     ceQsum/=norm;
1030     fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
1031     fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
1032   } else {
1033     ceQmax=0;
1034     ceTime=0;
1035     ceRMS =0;
1036     ceQsum=0;
1037   }
1038   param[0] = ceQmax;
1039   param[1] = ceTime;
1040   param[2] = ceRMS;
1041   qSum     = ceQsum;
1042 }
1043 //_____________________________________________________________________
1044 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
1045 {
1046   //
1047   // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
1048   // and 'tplus' timebins after 'pos'
1049   //
1050   if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
1051   for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
1052     if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
1053   for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
1054     if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
1055       ++iTime2;
1056     if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
1057   }
1058   return kTRUE;
1059 }
1060 //_____________________________________________________________________
1061 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
1062 {
1063   //
1064     // Find local maxima on the pad signal and Histogram them
1065   //
1066   Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal
1067   Int_t   count       = 0;
1068
1069   for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
1070     if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
1071     if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
1072       if (count<maxima.GetNrows()){
1073         maxima.GetMatrixArray()[count++]=i;
1074         GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
1075         i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin  fPeakDetMinus+fPeakDetPlus-1
1076       }
1077     }
1078   }
1079 }
1080 //_____________________________________________________________________
1081 void AliTPCCalibCE::ProcessPad()
1082 {
1083   //
1084   //  Process data of current pad
1085   //
1086   FindPedestal();
1087
1088   TVectorF maxima(15);     // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
1089                              // + central electrode and possibly post peaks from the CE signal
1090                              // however if we are on a high noise pad a lot more peaks due to the noise might occur
1091   FindLocalMaxima(maxima);
1092   if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
1093
1094   UpdateCETimeRef();                       // update the time refenrence for the current sector
1095   if ( fCurrentCETimeRef<1e-30 ) return;      //return if we don't have time 0 info, eg if only one side has laser
1096   TVectorD param(3);
1097   Float_t  qSum;
1098   FindCESignal(param, qSum, maxima);
1099
1100   Double_t meanT  = param[1];
1101   Double_t sigmaT = param[2];
1102
1103     //Fill Event T0 counter
1104   (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
1105
1106     //Fill Q histogram
1107   GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
1108
1109     //Fill RMS histogram
1110   GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
1111
1112
1113     //Fill debugging info
1114   if ( GetStreamLevel()>0 ){
1115     (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
1116     (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
1117     (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
1118   }
1119
1120   ResetPad();
1121 }
1122 //_____________________________________________________________________
1123 void AliTPCCalibCE::EndEvent()
1124 {
1125   //  Process data of current pad
1126   //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
1127   //  before the EndEvent function to set the event timestamp and number!!!
1128   //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
1129   //  function was called
1130   if (!fProcessOld) {
1131     if (fProcessNew) ++fNevents;
1132     return;
1133   }
1134
1135   //check if last pad has allready been processed, if not do so
1136   if ( fMaxTimeBin>-1 ) ProcessPad();
1137
1138   AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
1139
1140   TVectorD param(3);
1141   TMatrixD dummy(3,3);
1142 //    TVectorF vMeanTime(72);
1143 //    TVectorF vMeanQ(72);
1144   AliTPCCalROC *calIroc=new AliTPCCalROC(0);
1145   AliTPCCalROC *calOroc=new AliTPCCalROC(36);
1146
1147   //find mean time0 offset for side A and C
1148   //use only orocs due to the better statistics
1149   Double_t time0Side[2];       //time0 for side A:0 and C:1
1150   Double_t time0SideCount[2];  //time0 counter for side A:0 and C:1
1151   time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
1152   for ( Int_t iSec = 36; iSec<72; ++iSec ){
1153     time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
1154     time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
1155   }
1156   if ( time0SideCount[0] >0  )
1157     time0Side[0]/=time0SideCount[0];
1158   if ( time0SideCount[1] >0 )
1159     time0Side[1]/=time0SideCount[1];
1160     // end find time0 offset
1161   AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1162   Int_t nSecMeanT=0;
1163   //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1164   for ( Int_t iSec = 0; iSec<72; ++iSec ){
1165     AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1166     //find median and then calculate the mean around it
1167     TH1S *hMeanT    = GetHistoTmean(iSec); //histogram with local maxima position information
1168     if ( !hMeanT ) continue;
1169     //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1170     if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1171       hMeanT->Reset();
1172       AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1173       continue;
1174     }
1175
1176     Double_t entries = hMeanT->GetEffectiveEntries();
1177     Double_t sum     = 0;
1178     Short_t *arr     = hMeanT->GetArray()+1;
1179     Int_t ibin=0;
1180     for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1181       sum+=arr[ibin];
1182       if ( sum>=(entries/2.) ) break;
1183     }
1184     Int_t delta = 4;
1185     Int_t firstBin = fFirstTimeBin+ibin-delta;
1186     Int_t lastBin  = fFirstTimeBin+ibin+delta;
1187     if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1188     if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
1189     Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1190
1191         // check boundaries for ebye info of mean time
1192     TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1193     Int_t vSize=vMeanTime->GetNrows();
1194     if ( vSize < fNevents+1 ){
1195       vMeanTime->ResizeTo(vSize+100);
1196     }
1197
1198     // store mean time for the readout sides
1199     vSize=fVTime0SideA.GetNrows();
1200     if ( vSize < fNevents+1 ){
1201       fVTime0SideA.ResizeTo(vSize+100);
1202       fVTime0SideC.ResizeTo(vSize+100);
1203     }
1204     fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1205     fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
1206
1207     vMeanTime->GetMatrixArray()[fNevents]=median;
1208     nSecMeanT++;
1209     // end find median
1210
1211     TVectorF *vTimes = GetPadTimesEvent(iSec);
1212     if ( !vTimes ) continue;                     //continue if no time information for this sector is available
1213
1214     AliTPCCalROC calIrocOutliers(0);
1215     AliTPCCalROC calOrocOutliers(36);
1216
1217     // calculate mean Q of the sector
1218     TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1219     vSize=vMeanQ->GetNrows();
1220     if ( vSize < fNevents+1 ){
1221       vMeanQ->ResizeTo(vSize+100);
1222     }
1223     Float_t meanQ = 0;
1224     if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1225     vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1226
1227     for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1228       Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
1229
1230             //set values for temporary roc calibration class
1231       if ( iSec < 36 ) {
1232         calIroc->SetValue(iChannel, time);
1233         if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
1234
1235       } else {
1236         calOroc->SetValue(iChannel, time);
1237         if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
1238       }
1239
1240       if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1241         // test that we really found the CE signal reliably
1242         if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1243           GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
1244
1245
1246
1247             //-------------------------------  Debug start  ------------------------------
1248       if ( GetStreamLevel()>0 ){
1249         TTreeSRedirector *streamer=GetDebugStreamer();
1250         if (streamer){
1251           Int_t row=0;
1252           Int_t pad=0;
1253           Int_t padc=0;
1254
1255           Float_t q   = (*GetPadQEvent(iSec))[iChannel];
1256           Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1257
1258           UInt_t channel=iChannel;
1259           Int_t sector=iSec;
1260
1261           while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1262           pad = channel-fROC->GetRowIndexes(sector)[row];
1263           padc = pad-(fROC->GetNPads(sector,row)/2);
1264
1265 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1266 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
1267 //                                  fLastTimeBin-fFirstTimeBin,
1268 //                                  fFirstTimeBin,fLastTimeBin);
1269 //              h1->SetDirectory(0);
1270         //
1271 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1272 //                  h1->Fill(i,fPadSignal(i));
1273
1274           Double_t t0Sec = 0;
1275           if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1276             t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1277           Double_t t0Side = time0Side[(iSec/18)%2];
1278           (*streamer) << "DataPad" <<
1279             "Event=" << fNevents <<
1280             "RunNumber=" << fRunNumber <<
1281             "TimeStamp="   << fTimeStamp <<
1282             "Sector="<< sector <<
1283             "Row="   << row<<
1284             "Pad="   << pad <<
1285             "PadC="  << padc <<
1286             "PadSec="<< channel <<
1287             "Time0Sec="  << t0Sec <<
1288             "Time0Side=" << t0Side <<
1289             "Time="  << time <<
1290             "RMS="   << rms <<
1291             "Sum="   << q <<
1292             "MeanQ=" << meanQ <<
1293         //                  "hist.=" << h1 <<
1294             "\n";
1295
1296     //          delete h1;
1297         }
1298       }
1299       //-----------------------------  Debug end  ------------------------------
1300     }// end channel loop
1301
1302
1303     //do fitting now only in debug mode
1304     if (GetDebugLevel()>0){
1305       TVectorD paramPol1(3);
1306       TVectorD paramPol2(6);
1307       TMatrixD matPol1(3,3);
1308       TMatrixD matPol2(6,6);
1309       Float_t  chi2Pol1=0;
1310       Float_t  chi2Pol2=0;
1311
1312       if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1313         if ( iSec < 36 ){
1314           calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1315           calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1316         } else {
1317           calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1318           calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1319         }
1320
1321         GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1322         GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1323       }
1324
1325   //-------------------------------  Debug start  ------------------------------
1326       if ( GetStreamLevel()>0 ){
1327         TTreeSRedirector *streamer=GetDebugStreamer();
1328         if ( streamer ) {
1329           (*streamer) << "DataRoc" <<
1330 //    "Event=" << fEvent <<
1331             "RunNumber=" << fRunNumber <<
1332             "TimeStamp="   << fTimeStamp <<
1333             "Sector="<< iSec <<
1334             "hMeanT.=" << hMeanT <<
1335             "median=" << median <<
1336             "paramPol1.=" << &paramPol1 <<
1337             "paramPol2.=" << &paramPol2 <<
1338             "matPol1.="   << &matPol1 <<
1339             "matPol2.="   << &matPol2 <<
1340             "chi2Pol1="   << chi2Pol1 <<
1341             "chi2Pol2="   << chi2Pol2 <<
1342             "\n";
1343         }
1344       }
1345     }
1346         //-------------------------------  Debug end  ------------------------------
1347     hMeanT->Reset();
1348   }// end sector loop
1349     //return if no sector has a valid mean time
1350   if ( nSecMeanT == 0 ) return;
1351
1352
1353 //    fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1354 //    fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1355   if ( fVEventTime.GetNrows() < fNevents+1 ) {
1356     fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1357     fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1358   }
1359   fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1360   fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1361
1362   ++fNevents;
1363   fOldRunNumber = fRunNumber;
1364
1365   delete calIroc;
1366   delete calOroc;
1367   AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1368 }
1369 //_____________________________________________________________________
1370 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1371                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
1372                                   const Char_t *type, Bool_t force)
1373 {
1374     //
1375     // return pointer to TH2S histogram of 'type'
1376     // if force is true create a new histogram if it doesn't exist allready
1377     //
1378     if ( !force || arr->UncheckedAt(sector) )
1379       return (TH2S*)arr->UncheckedAt(sector);
1380
1381     // if we are forced and histogram doesn't exist yet create it
1382     // new histogram with Q calib information. One value for each pad!
1383     TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1384                           nbinsY, ymin, ymax,
1385                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1386     hist->SetDirectory(0);
1387     arr->AddAt(hist,sector);
1388     return hist;
1389 }
1390 //_____________________________________________________________________
1391 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1392 {
1393     //
1394     // return pointer to T0 histogram
1395     // if force is true create a new histogram if it doesn't exist allready
1396     //
1397     TObjArray *arr = &fHistoT0Array;
1398     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1399 }
1400 //_____________________________________________________________________
1401 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1402 {
1403     //
1404     // return pointer to Q histogram
1405     // if force is true create a new histogram if it doesn't exist allready
1406     //
1407     TObjArray *arr = &fHistoQArray;
1408     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1409 }
1410 //_____________________________________________________________________
1411 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1412 {
1413     //
1414     // return pointer to Q histogram
1415     // if force is true create a new histogram if it doesn't exist allready
1416     //
1417     TObjArray *arr = &fHistoRMSArray;
1418     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1419 }
1420 //_____________________________________________________________________
1421 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1422                               const Char_t *type, Bool_t force)
1423 {
1424     //
1425     // return pointer to TH1S histogram
1426     // if force is true create a new histogram if it doesn't exist allready
1427     //
1428     if ( !force || arr->UncheckedAt(sector) )
1429       return (TH1S*)arr->UncheckedAt(sector);
1430
1431     // if we are forced and histogram doesn't yes exist create it
1432     // new histogram with calib information. One value for each pad!
1433     TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1434                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1435     hist->SetDirectory(0);
1436     arr->AddAt(hist,sector);
1437     return hist;
1438 }
1439 //_____________________________________________________________________
1440 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1441 {
1442     //
1443     // return pointer to Q histogram
1444     // if force is true create a new histogram if it doesn't exist allready
1445     //
1446     TObjArray *arr = &fHistoTmean;
1447     return GetHisto(sector, arr, "LastTmean", force);
1448 }
1449 //_____________________________________________________________________
1450 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1451 {
1452     //
1453     // return pointer to Pad Info from 'arr' for the current event and sector
1454     // if force is true create it if it doesn't exist allready
1455     //
1456     if ( !force || arr->UncheckedAt(sector) )
1457         return (TVectorF*)arr->UncheckedAt(sector);
1458
1459     TVectorF *vect = new TVectorF(size);
1460     arr->AddAt(vect,sector);
1461     return vect;
1462 }
1463 //_____________________________________________________________________
1464 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1465 {
1466     //
1467     // return pointer to Pad Times Array for the current event and sector
1468     // if force is true create it if it doesn't exist allready
1469     //
1470     TObjArray *arr = &fPadTimesArrayEvent;
1471     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1472 }
1473 //_____________________________________________________________________
1474 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1475 {
1476     //
1477     // return pointer to Pad Q Array for the current event and sector
1478     // if force is true create it if it doesn't exist allready
1479     // for debugging purposes only
1480     //
1481
1482     TObjArray *arr = &fPadQArrayEvent;
1483     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1484 }
1485 //_____________________________________________________________________
1486 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1487 {
1488     //
1489     // return pointer to Pad RMS Array for the current event and sector
1490     // if force is true create it if it doesn't exist allready
1491     // for debugging purposes only
1492     //
1493     TObjArray *arr = &fPadRMSArrayEvent;
1494     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1495 }
1496 //_____________________________________________________________________
1497 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1498 {
1499     //
1500     // return pointer to Pad RMS Array for the current event and sector
1501     // if force is true create it if it doesn't exist allready
1502     // for debugging purposes only
1503     //
1504     TObjArray *arr = &fPadPedestalArrayEvent;
1505     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1506 }
1507 //_____________________________________________________________________
1508 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1509 {
1510     //
1511     // return pointer to the EbyE info of the mean arrival time for 'sector'
1512     // if force is true create it if it doesn't exist allready
1513     //
1514     TObjArray *arr = &fTMeanArrayEvent;
1515     return GetVectSector(sector,arr,100,force);
1516 }
1517 //_____________________________________________________________________
1518 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1519 {
1520     //
1521     // return pointer to the EbyE info of the mean arrival time for 'sector'
1522     // if force is true create it if it doesn't exist allready
1523     //
1524     TObjArray *arr = &fQMeanArrayEvent;
1525     return GetVectSector(sector,arr,100,force);
1526 }
1527 //_____________________________________________________________________
1528 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1529 {
1530     //
1531     // return pointer to ROC Calibration
1532     // if force is true create a new histogram if it doesn't exist allready
1533     //
1534     if ( !force || arr->UncheckedAt(sector) )
1535         return (AliTPCCalROC*)arr->UncheckedAt(sector);
1536
1537     // if we are forced and histogram doesn't yes exist create it
1538
1539     // new AliTPCCalROC for T0 information. One value for each pad!
1540     AliTPCCalROC *croc = new AliTPCCalROC(sector);
1541     arr->AddAt(croc,sector);
1542     return croc;
1543 }
1544 //_____________________________________________________________________
1545 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1546 {
1547     //
1548     // return pointer to Time 0 ROC Calibration
1549     // if force is true create a new histogram if it doesn't exist allready
1550     //
1551     TObjArray *arr = &fCalRocArrayT0;
1552     return GetCalRoc(sector, arr, force);
1553 }
1554 //_____________________________________________________________________
1555 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1556 {
1557     //
1558     // return pointer to the error of Time 0 ROC Calibration
1559     // if force is true create a new histogram if it doesn't exist allready
1560     //
1561     TObjArray *arr = &fCalRocArrayT0Err;
1562     return GetCalRoc(sector, arr, force);
1563 }
1564 //_____________________________________________________________________
1565 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1566 {
1567     //
1568     // return pointer to T0 ROC Calibration
1569     // if force is true create a new histogram if it doesn't exist allready
1570     //
1571     TObjArray *arr = &fCalRocArrayQ;
1572     return GetCalRoc(sector, arr, force);
1573 }
1574 //_____________________________________________________________________
1575 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1576 {
1577     //
1578     // return pointer to signal width ROC Calibration
1579     // if force is true create a new histogram if it doesn't exist allready
1580     //
1581     TObjArray *arr = &fCalRocArrayRMS;
1582     return GetCalRoc(sector, arr, force);
1583 }
1584 //_____________________________________________________________________
1585 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1586 {
1587     //
1588     // return pointer to Outliers
1589     // if force is true create a new histogram if it doesn't exist allready
1590     //
1591     TObjArray *arr = &fCalRocArrayOutliers;
1592     return GetCalRoc(sector, arr, force);
1593 }
1594 //_____________________________________________________________________
1595 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1596 {
1597     //
1598     // return pointer to TObjArray of fit parameters
1599     // if force is true create a new histogram if it doesn't exist allready
1600     //
1601     if ( !force || arr->UncheckedAt(sector) )
1602         return (TObjArray*)arr->UncheckedAt(sector);
1603
1604     // if we are forced and array doesn't yes exist create it
1605
1606     // new TObjArray for parameters
1607     TObjArray *newArr = new TObjArray;
1608     arr->AddAt(newArr,sector);
1609     return newArr;
1610 }
1611 //_____________________________________________________________________
1612 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1613 {
1614     //
1615     // return pointer to TObjArray of fit parameters from plane fit
1616     // if force is true create a new histogram if it doesn't exist allready
1617     //
1618     TObjArray *arr = &fParamArrayEventPol1;
1619     return GetParamArray(sector, arr, force);
1620 }
1621 //_____________________________________________________________________
1622 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1623 {
1624     //
1625     // return pointer to TObjArray of fit parameters from parabola fit
1626     // if force is true create a new histogram if it doesn't exist allready
1627     //
1628     TObjArray *arr = &fParamArrayEventPol2;
1629     return GetParamArray(sector, arr, force);
1630 }
1631
1632 //_____________________________________________________________________
1633 void AliTPCCalibCE::CreateDVhist()
1634 {
1635   //
1636   // Setup the THnSparse for the drift velocity determination
1637   //
1638
1639   //HnSparse bins
1640   //roc, row, pad, timebin, timestamp
1641   TTimeStamp begin(2010,01,01,0,0,0);
1642   TTimeStamp end(2035,01,01,0,0,0);
1643   Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1644
1645   Int_t    bins[kHnBinsDV] = { 72,  96,  140,  1030, nbinsTime};
1646   Double_t xmin[kHnBinsDV] = { 0.,  0.,   0.,    0., (Double_t)begin.GetSec()};
1647   Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1648
1649   fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1650   fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1651   fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1652   fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1653   fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1654   fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
1655
1656 }
1657
1658 //_____________________________________________________________________
1659 void AliTPCCalibCE::ResetEvent()
1660 {
1661     //
1662     //  Reset global counters  -- Should be called before each event is processed
1663     //
1664     fLastSector=-1;
1665     fCurrentSector=-1;
1666     fCurrentRow=-1;
1667     fCurrentChannel=-1;
1668
1669     ResetPad();
1670
1671     fPadTimesArrayEvent.Delete();
1672     fPadQArrayEvent.Delete();
1673     fPadRMSArrayEvent.Delete();
1674     fPadPedestalArrayEvent.Delete();
1675
1676     for ( Int_t i=0; i<72; ++i ){
1677         fVTime0Offset.GetMatrixArray()[i]=0;
1678         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1679         fVMeanQ.GetMatrixArray()[i]=0;
1680         fVMeanQCounter.GetMatrixArray()[i]=0;
1681     }
1682 }
1683 //_____________________________________________________________________
1684 void AliTPCCalibCE::ResetPad()
1685 {
1686     //
1687     //  Reset pad infos -- Should be called after a pad has been processed
1688     //
1689     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1690       fPadSignal[i] = 0;
1691     fMaxTimeBin   = -1;
1692     fMaxPadSignal = -1;
1693     fPadPedestal  = -1;
1694     fPadNoise     = -1;
1695 }
1696 //_____________________________________________________________________
1697 void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
1698 {
1699   //
1700   //  Merge ce to the current AliTPCCalibCE
1701   //
1702   MergeBase(ce);
1703   Int_t nCEevents = ce->GetNeventsProcessed();
1704
1705   if (fProcessOld&&ce->fProcessOld){
1706   //merge histograms
1707     for (Int_t iSec=0; iSec<72; ++iSec){
1708       TH2S *hRefQmerge   = ce->GetHistoQ(iSec);
1709       TH2S *hRefT0merge  = ce->GetHistoT0(iSec);
1710       TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1711
1712
1713       if ( hRefQmerge ){
1714         TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1715         TH2S *hRefQ   = GetHistoQ(iSec);
1716         if ( hRefQ ) hRefQ->Add(hRefQmerge);
1717         else {
1718           TH2S *hist = new TH2S(*hRefQmerge);
1719           hist->SetDirectory(0);
1720           fHistoQArray.AddAt(hist, iSec);
1721         }
1722         hRefQmerge->SetDirectory(dir);
1723       }
1724       if ( hRefT0merge ){
1725         TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1726         TH2S *hRefT0  = GetHistoT0(iSec);
1727         if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1728         else {
1729           TH2S *hist = new TH2S(*hRefT0merge);
1730           hist->SetDirectory(0);
1731           fHistoT0Array.AddAt(hist, iSec);
1732         }
1733         hRefT0merge->SetDirectory(dir);
1734       }
1735       if ( hRefRMSmerge ){
1736         TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1737         TH2S *hRefRMS = GetHistoRMS(iSec);
1738         if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1739         else {
1740           TH2S *hist = new TH2S(*hRefRMSmerge);
1741           hist->SetDirectory(0);
1742           fHistoRMSArray.AddAt(hist, iSec);
1743         }
1744         hRefRMSmerge->SetDirectory(dir);
1745       }
1746
1747     }
1748
1749     // merge time information
1750
1751
1752     for (Int_t iSec=0; iSec<72; ++iSec){
1753       TObjArray *arrPol1CE  = ce->GetParamArrayPol1(iSec);
1754       TObjArray *arrPol2CE  = ce->GetParamArrayPol2(iSec);
1755       TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1756       TVectorF *vMeanQCE    = ce->GetQMeanEvents(iSec);
1757
1758       TObjArray *arrPol1  = 0x0;
1759       TObjArray *arrPol2  = 0x0;
1760       TVectorF *vMeanTime = 0x0;
1761       TVectorF *vMeanQ    = 0x0;
1762
1763   //resize arrays
1764       if ( arrPol1CE && arrPol2CE ){
1765         arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1766         arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1767         arrPol1->Expand(fNevents+nCEevents);
1768         arrPol2->Expand(fNevents+nCEevents);
1769       }
1770       if ( vMeanTimeCE && vMeanQCE ){
1771         vMeanTime = GetTMeanEvents(iSec,kTRUE);
1772         vMeanQ    = GetQMeanEvents(iSec,kTRUE);
1773         vMeanTime->ResizeTo(fNevents+nCEevents);
1774         vMeanQ->ResizeTo(fNevents+nCEevents);
1775       }
1776
1777       for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1778         if ( arrPol1CE && arrPol2CE ){
1779           TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1780           TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1781           if ( paramPol1 && paramPol2 ){
1782             GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1783             GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1784           }
1785         }
1786         if ( vMeanTimeCE && vMeanQCE ){
1787           vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1788           vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1789         }
1790       }
1791     }
1792
1793
1794
1795     const TVectorD&  eventTimes  = ce->fVEventTime;
1796     const TVectorD&  eventIds    = ce->fVEventNumber;
1797     const TVectorF&  time0SideA  = ce->fVTime0SideA;
1798     const TVectorF&  time0SideC  = ce->fVTime0SideC;
1799     fVEventTime.ResizeTo(fNevents+nCEevents);
1800     fVEventNumber.ResizeTo(fNevents+nCEevents);
1801     fVTime0SideA.ResizeTo(fNevents+nCEevents);
1802     fVTime0SideC.ResizeTo(fNevents+nCEevents);
1803
1804     for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1805       Double_t evTime     = eventTimes.GetMatrixArray()[iEvent];
1806       Double_t evId       = eventIds.GetMatrixArray()[iEvent];
1807       Float_t  t0SideA    = time0SideA.GetMatrixArray()[iEvent];
1808       Float_t  t0SideC    = time0SideC.GetMatrixArray()[iEvent];
1809
1810       fVEventTime.GetMatrixArray()[fNevents+iEvent]   = evTime;
1811       fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1812       fVTime0SideA.GetMatrixArray()[fNevents+iEvent]  = t0SideA;
1813       fVTime0SideC.GetMatrixArray()[fNevents+iEvent]  = t0SideC;
1814     }
1815   }
1816
1817   if (fProcessNew&&ce->fProcessNew) {
1818     if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1819       AliError("Number of bursts in the instances to merge are different. No merging done!");
1820     } else {
1821       for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1822         THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1823         THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1824         if (h && hce) h->Add(hce);
1825         else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
1826       }
1827       //TODO: What to do with fTimeBursts???
1828     }
1829   }
1830
1831   fNevents+=nCEevents; //increase event counter
1832 }
1833
1834 //_____________________________________________________________________
1835 Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1836 {
1837   //
1838   // Merge all objects of this type in list
1839   //
1840
1841   Long64_t nmerged=1;
1842
1843   TIter next(list);
1844   AliTPCCalibCE *ce=0;
1845   TObject *o=0;
1846
1847   while ( (o=next()) ){
1848     ce=dynamic_cast<AliTPCCalibCE*>(o);
1849     if (ce){
1850       Merge(ce);
1851       ++nmerged;
1852     }
1853   }
1854
1855   return nmerged;
1856 }
1857
1858 //_____________________________________________________________________
1859 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1860 {
1861   //
1862   // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1863   // or side (-1: A-Side, -2: C-Side)
1864   // xVariable:    0-event time, 1-event id, 2-internal event counter
1865   // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1866   // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1867   //                               0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1868   //                               not used for mean time and mean Q )
1869   // for an example see class description at the beginning
1870   //
1871
1872   TVectorD *xVar = 0x0;
1873   TObjArray *aType = 0x0;
1874   Int_t npoints=0;
1875
1876   // sanity checks
1877   if ( (sector<-2)   || (sector>71)   ) return 0x0;  //sector outside valid range
1878   if ( (xVariable<0) || (xVariable>2) ) return 0x0;  //invalid x-variable
1879   if ( (fitType<0)   || (fitType>3)   ) return 0x0;  //invalid fit type
1880   if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1881   if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1882   if ( sector<0 && fitType!=2) return 0x0;  //for side wise information only mean time is available
1883
1884   if (sector>=0){
1885     if ( fitType==0 ){
1886       if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1887       aType = &fParamArrayEventPol1;
1888       if ( aType->At(sector)==0x0 ) return 0x0;
1889     }
1890     else if ( fitType==1 ){
1891       if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1892       aType = &fParamArrayEventPol2;
1893       if ( aType->At(sector)==0x0 ) return 0x0;
1894     }
1895
1896   }
1897   if ( xVariable == 0 ) xVar = &fVEventTime;
1898   if ( xVariable == 1 ) xVar = &fVEventNumber;
1899   if ( xVariable == 2 ) {
1900     xVar = new TVectorD(fNevents);
1901     for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1902   }
1903
1904   Double_t *x = new Double_t[fNevents];
1905   Double_t *y = new Double_t[fNevents];
1906
1907   for (Int_t ievent =0; ievent<fNevents; ++ievent){
1908     if ( fitType<2 ){
1909       TObjArray *events = (TObjArray*)(aType->At(sector));
1910       if ( events->GetSize()<=ievent ) break;
1911       TVectorD *v = (TVectorD*)(events->At(ievent));
1912       if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1913     } else if (fitType == 2) {
1914       Double_t xValue=(*xVar)[ievent];
1915       Double_t yValue=0;
1916       if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1917       else if (sector==-1) yValue=fVTime0SideA(ievent);
1918       else if (sector==-2) yValue=fVTime0SideC(ievent);
1919       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1920     }else if (fitType == 3) {
1921       Double_t xValue=(*xVar)[ievent];
1922       Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1923       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1924     }
1925   }
1926
1927   TGraph *gr = new TGraph(npoints);
1928     //sort xVariable increasing
1929   Int_t    *sortIndex = new Int_t[npoints];
1930   TMath::Sort(npoints,x,sortIndex, kFALSE);
1931   for (Int_t i=0;i<npoints;++i){
1932     gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1933   }
1934
1935
1936   if ( xVariable == 2 ) delete xVar;
1937   delete [] x;
1938   delete [] y;
1939   delete [] sortIndex;
1940   return gr;
1941 }
1942 //_____________________________________________________________________
1943 void AliTPCCalibCE::Analyse()
1944 {
1945   //
1946   //  Calculate calibration constants
1947   //
1948
1949   if (fProcessOld){
1950     TVectorD paramQ(3);
1951     TVectorD paramT0(3);
1952     TVectorD paramRMS(3);
1953     TMatrixD dummy(3,3);
1954
1955     Float_t channelCounter=0;
1956     fMeanT0rms=0;
1957     fMeanQrms=0;
1958     fMeanRMSrms=0;
1959
1960     for (Int_t iSec=0; iSec<72; ++iSec){
1961       TH2S *hT0 = GetHistoT0(iSec);
1962       if (!hT0 ) continue;
1963
1964       AliTPCCalROC *rocQ     = GetCalRocQ  (iSec,kTRUE);
1965       AliTPCCalROC *rocT0    = GetCalRocT0 (iSec,kTRUE);
1966       AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1967       AliTPCCalROC *rocRMS   = GetCalRocRMS(iSec,kTRUE);
1968       AliTPCCalROC *rocOut   = GetCalRocOutliers(iSec,kTRUE);
1969
1970       TH2S *hQ   = GetHistoQ(iSec);
1971       TH2S *hRMS = GetHistoRMS(iSec);
1972
1973       Short_t *arrayhQ   = hQ->GetArray();
1974       Short_t *arrayhT0  = hT0->GetArray();
1975       Short_t *arrayhRMS = hRMS->GetArray();
1976
1977       UInt_t nChannels = fROC->GetNChannels(iSec);
1978
1979   //debug
1980       Int_t row=0;
1981       Int_t pad=0;
1982       Int_t padc=0;
1983   //! debug
1984
1985       for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1986
1987
1988         Float_t cogTime0 = -1000;
1989         Float_t cogQ     = -1000;
1990         Float_t cogRMS   = -1000;
1991         Float_t cogOut   = 0;
1992         Float_t rms      = 0;
1993         Float_t rmsT0    = 0;
1994
1995
1996         Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1997         Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1998         Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1999
2000         cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2001         fMeanQrms+=rms;
2002         cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2003         fMeanT0rms+=rmsT0;
2004         cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2005         fMeanRMSrms+=rms;
2006         channelCounter++;
2007
2008       /*
2009              //outlier specifications
2010          if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2011          cogOut = 1;
2012           cogTime0 = 0;
2013           cogQ     = 0;
2014           cogRMS   = 0;
2015       }
2016       */
2017         rocQ->SetValue(iChannel, cogQ*cogQ);
2018         rocT0->SetValue(iChannel, cogTime0);
2019         rocT0Err->SetValue(iChannel, rmsT0);
2020         rocRMS->SetValue(iChannel, cogRMS);
2021         rocOut->SetValue(iChannel, cogOut);
2022
2023
2024       //debug
2025         if ( GetStreamLevel() > 0 ){
2026           TTreeSRedirector *streamer=GetDebugStreamer();
2027           if ( streamer ) {
2028
2029             while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2030             pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2031             padc = pad-(fROC->GetNPads(iSec,row)/2);
2032
2033             (*streamer) << "DataEnd" <<
2034               "Sector="  << iSec      <<
2035               "Pad="     << pad       <<
2036               "PadC="    << padc      <<
2037               "Row="     << row       <<
2038               "PadSec="  << iChannel   <<
2039               "Q="       << cogQ      <<
2040               "T0="      << cogTime0  <<
2041               "RMS="     << cogRMS    <<
2042               "\n";
2043           }
2044         }
2045       //! debug
2046
2047       }
2048
2049     }
2050     if ( channelCounter>0 ){
2051       fMeanT0rms/=channelCounter;
2052       fMeanQrms/=channelCounter;
2053       fMeanRMSrms/=channelCounter;
2054     }
2055     //   if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2056     //    delete fDebugStreamer;
2057     //    fDebugStreamer = 0x0;
2058     fVEventTime.ResizeTo(fNevents);
2059     fVEventNumber.ResizeTo(fNevents);
2060     fVTime0SideA.ResizeTo(fNevents);
2061     fVTime0SideC.ResizeTo(fNevents);
2062   }
2063
2064   if (fProcessNew&&fAnalyseNew){
2065     AnalyseTrack();
2066     for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2067       THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2068       h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2069     }
2070   }
2071 }
2072
2073
2074
2075
2076 //
2077 // New functions that also use the laser tracks
2078 //
2079
2080
2081
2082 //_____________________________________________________________________
2083 void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2084 {
2085   //
2086   //Find the local maximums for the projections to each axis
2087   //
2088
2089   //find laser layer positoins
2090   fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2091   FindLaserLayers();
2092   THnSparse *hProj=fHnDrift;
2093   Double_t posCE[4]={0.,0.,0.,0.};
2094   Double_t widthCE[4]={0.,0.,0.,0.};
2095
2096 //   if(fPeaks[4]!=0){
2097   // find central electrode position once more, separately for IROC, OROC, A-, C-Side
2098
2099   for (Int_t i=0; i<4; ++i){
2100     hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2101     TH1 *h=fHnDrift->Projection(3);
2102     h->GetXaxis()->SetRangeUser(fPeaks[4]-fPeakWidths[4],fPeaks[4]+fPeakWidths[4]);
2103     Int_t nbinMax=h->GetMaximumBin();
2104     Double_t maximum=h->GetMaximum();
2105 //     Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2106 //     if (nbinMax<700||maximum<maxExpected) continue;
2107     Double_t xbinMax=h->GetBinCenter(nbinMax);
2108     TF1 fgaus("gaus","gaus",xbinMax-10,xbinMax+10);
2109     fgaus.SetParameters(maximum,xbinMax,2);
2110     fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2111     fgaus.SetParLimits(2,0.2,4.);
2112     h->Fit(&fgaus,"RQN");
2113 //     Double_t deltaX=4*fgaus.GetParameter(2);
2114 //     xbinMax=fgaus.GetParameter(1);
2115     delete h;
2116     posCE[i]=fgaus.GetParameter(1);
2117     widthCE[i]=4*fgaus.GetParameter(2);
2118     hProj->GetAxis(0)->SetRangeUser(0,72);
2119   }
2120 //   }
2121   //Current drift velocity
2122   Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2123 //   cout<<"vdrift="<<vdrift<<endl;
2124
2125   AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2126   //loop over all entries in the histogram
2127   Int_t coord[5];
2128   for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2129     //get entry position and content
2130     Double_t adc=hProj->GetBinContent(ichunk,coord);
2131
2132
2133     Int_t sector = coord[0]-1;
2134     Int_t row    = coord[1]-1;
2135     Int_t pad    = coord[2]-1;
2136     Int_t timeBin= coord[3]-1;
2137     Double_t time   = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2138     Int_t side   = (sector/18)%2;
2139 //     return;
2140 //     fPeaks[4]=(UInt_t)posCE[sector/18];
2141 //     fPeakWidths[4]=(UInt_t)widthCE[sector/18];
2142
2143     //cuts
2144     if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2145     if (adc < 5 ) continue;
2146     if (IsEdgePad(sector,row,pad)) continue;
2147 //     if (!IsPeakInRange(timeBin)) continue;
2148 //     if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2149 //     if (isCE&&(sector!=0)) continue;
2150
2151     Int_t padmin=-2, padmax=2;
2152     Int_t timemin=-2, timemax=2;
2153     Int_t minsumperpad=2;
2154     //CE or laser tracks
2155     Bool_t isCE=kFALSE;
2156     if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2157       isCE=kTRUE;
2158       padmin=0;
2159       padmax=0;
2160       timemin=-3;
2161       timemax=7;
2162     }
2163
2164     //
2165     // Find local maximum and cogs
2166     //
2167     Bool_t isMaximum=kTRUE;
2168     Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2169     Double_t cogY=0, rmsY=0;
2170     Int_t npart=0;
2171
2172     // for position calculation use
2173     for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2174       Float_t lxyz[3];
2175       fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2176
2177       for(Int_t itime=timemin;itime<=timemax;++itime){
2178
2179         Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2180         Double_t val=hProj->GetBinContent(a);
2181         totalmass+=val;
2182
2183         tbcm +=(timeBin+itime)*val;
2184         padcm+=(pad+ipad)*val;
2185         cogY +=lxyz[1]*val;
2186
2187         rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2188         rmspad+=(pad+ipad)*(pad+ipad)*val;
2189         rmsY  +=lxyz[1]*lxyz[1]*val;
2190
2191         if (val>0) ++npart;
2192         if (val>adc) {
2193           isMaximum=kFALSE;
2194           break;
2195         }
2196       }
2197       if (!isMaximum)  break;
2198     }
2199
2200     if (!isMaximum||!npart)  continue;
2201     if (totalmass<npart*minsumperpad) continue;
2202     if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2203                                     //for CE we want only one pad by construction
2204
2205     tbcm/=totalmass;
2206     padcm/=totalmass;
2207     cogY/=totalmass;
2208
2209     rmstb/=totalmass;
2210     rmspad/=totalmass;
2211     rmsY/=totalmass;
2212
2213     rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2214     rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2215     rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2216
2217     Int_t cog=TMath::Nint(padcm);
2218
2219     // timebin --> z position
2220     Float_t zlength=fParam->GetZLength(side);
2221 //     Float_t timePos=tbcm+(1000-fPeaks[4])
2222     // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2223     Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2224
2225     // local to global transformation--> x and y positions
2226     Float_t padlxyz[3];
2227     fROC->GetPositionLocal(sector,row,pad,padlxyz);
2228
2229     Double_t gxyz[3]={padlxyz[0],cogY,gz};
2230     Double_t lxyz[3]={padlxyz[0],cogY,gz};
2231     Double_t igxyz[3]={0,0,0};
2232     AliTPCTransform t1;
2233     t1.RotatedGlobal2Global(sector,gxyz);
2234
2235     Double_t mindist=0;
2236     Int_t trackID=-1;
2237     Int_t trackID2=-1;
2238
2239     //find track id and index of the position in the track (row)
2240     Int_t index=0;
2241     if (!isCE){
2242       index=row+(sector>35)*fROC->GetNRows(0);
2243       trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2244     } else {
2245       trackID=336+((sector/18)%2);
2246       index= fROC->GetRowIndexes(sector)[row]+pad; //  global pad position in sector
2247       if (sector<36) {
2248         index+=(sector%18)*fROC->GetNChannels(sector);
2249       } else {
2250         index+=18*fROC->GetNChannels(0);
2251         index+=(sector%18)*fROC->GetNChannels(sector);
2252       }
2253       //TODO: find out about the multiple peaks in the CE
2254 //       mindist=TMath::Abs(fPeaks[4]-tbcm);
2255       mindist=1.;
2256     }
2257
2258     // fill track vectors
2259     if (trackID>0){
2260       AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2261       Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2262
2263
2264 //      travel time effect of light includes
2265
2266       Double_t raylength=ltr->GetRayLength();
2267       Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2268       ltr->GetXYZ(globmir);
2269       if(trackID<336){
2270         if(side==0){
2271           gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2272                                        +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2273                                        +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2274         }
2275         else {
2276           gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2277                                        +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2278                                        +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2279         }
2280       }
2281
2282       if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2283         ltr->fVecSec->GetMatrixArray()[index]=sector;
2284         ltr->fVecP2->GetMatrixArray()[index]=0;
2285         ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2286         ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2287         ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2288         ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2289         ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2290         ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2291         ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2292 //         ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2293       }
2294       TObjArray *arr=AliTPCLaserTrack::GetTracks();
2295       ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2296       igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2297       igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2298       igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2299     }
2300
2301
2302     if (fStreamLevel>4){
2303       (*GetDebugStreamer()) << "clusters" <<
2304         "run="   << fRunNumber <<
2305         "timestamp=" << timestamp <<
2306         "burst="     << burst     <<
2307         "side="      << side      <<
2308         "sec="       << sector    <<
2309         "row="       << row       <<
2310         "pad="       << pad       <<
2311         "padCog="    << cog       <<
2312         "timebin="   << timeBin   <<
2313         "cogCE="     << posCE[sector/18] <<
2314         "withCE="    << widthCE[sector/18] <<
2315         "index="     << index     <<
2316
2317         "padcm="     << padcm     <<
2318         "rmspad="    << rmspad    <<
2319
2320         "cogtb="     << tbcm      <<
2321         "rmstb="     << rmstb     <<
2322
2323         "npad="      << npart     <<
2324
2325         "lx="        << padlxyz[0]<<
2326         "ly="        << cogY      <<
2327         "lypad="     << padlxyz[1]<<
2328         "rmsY="      << rmsY      <<
2329
2330         "gx="        << gxyz[0]   <<
2331         "gy="        << gxyz[1]   <<
2332         "gz="        << gxyz[2]   <<
2333
2334         "igx="        << igxyz[0] <<
2335         "igy="        << igxyz[1] <<
2336         "igz="        << igxyz[2] <<
2337
2338         "mind="      << mindist   <<
2339         "max="       << adc       <<
2340         "trackid="   << trackID   <<
2341         "trackid2="   << trackID2   <<
2342         "npart="     << npart     <<
2343         "\n";
2344     } // end stream levelmgz.fElements
2345
2346   }
2347
2348 }
2349
2350 //_____________________________________________________________________
2351 void AliTPCCalibCE::AnalyseTrack()
2352 {
2353   //
2354   //  Analyse the tracks
2355   //
2356
2357
2358   AliTPCLaserTrack::LoadTracks();
2359 //   AliTPCParam *param=0x0;
2360 //   //cdb run number
2361 //   AliCDBManager *man=AliCDBManager::Instance();
2362 //   if (man->GetDefaultStorage()){
2363 //     AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2364 //     if (entry){
2365 //       entry->SetOwner(kTRUE);
2366 //       param = (AliTPCParam*)(entry->GetObject()->Clone());
2367 //     }
2368 //   }
2369 //   if (param){
2370 //     if (fParam) delete fParam;
2371 //     fParam=param;
2372 //   } else {
2373 //     AliError("Could not get updated AliTPCParam from OCDB!!!");
2374 //   }
2375
2376   //Measured and ideal laser tracks
2377   TObjArray* arrMeasured = SetupMeasured();
2378   TObjArray* arrIdeal    = AliTPCLaserTrack::GetTracks();
2379   AddCEtoIdeal(arrIdeal);
2380
2381   //find bursts and loop over them
2382   for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2383     Double_t timestamp=fTimeBursts[iburst];
2384     AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2385     fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2386     if (!fHnDrift) continue;
2387     UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2388     if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2389     fBinsLastAna[iburst]=entries;
2390
2391     for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2392 //     if (iburst==0) FindLaserLayers();
2393
2394     //reset laser tracks
2395     ResetMeasured(arrMeasured);
2396
2397     //find clusters and associate to the tracks
2398     FindLocalMaxima(arrMeasured, timestamp, iburst);
2399
2400     //calculate drift velocity
2401     CalculateDV(arrIdeal,arrMeasured,iburst);
2402
2403     //Dump information to file if requested
2404     if (fStreamLevel>2){
2405       printf("make tree\n");
2406       //laser track information
2407
2408       for (Int_t itrack=0; itrack<338; ++itrack){
2409         TObject *iltr=arrIdeal->UncheckedAt(itrack);
2410         TObject *mltr=arrMeasured->UncheckedAt(itrack);
2411         (*GetDebugStreamer()) << "tracks" <<
2412           "run="   << fRunNumber <<
2413           "time=" << timestamp <<
2414           "burst="<< iburst <<
2415           "iltr.=" << iltr <<
2416           "mltr.=" << mltr <<
2417           "\n";
2418       }
2419     }
2420   }
2421   if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2422 }
2423
2424 //_____________________________________________________________________
2425 Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2426                                       const Double_t *peakposloc, Int_t &itrackMin2)
2427 {
2428   //
2429   //  Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2430   //
2431
2432
2433   TObjArray *arr=AliTPCLaserTrack::GetTracks();
2434   TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2435   TVector3 vDir;
2436   TVector3 vSt;
2437
2438   Int_t firstbeam=0;
2439   Int_t lastbeam=336/2;
2440   if ( (sector/18)%2 ) {
2441     firstbeam=336/2;
2442     lastbeam=336;
2443   }
2444
2445   mindist=1000000;
2446   Int_t itrackMin=-1;
2447   for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2448     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);  //get the track
2449 //     if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2450     vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2451     Double_t deltaZ=ltr->GetZ()-peakpos[2];
2452     if (TMath::Abs(deltaZ)>40) continue;
2453     vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2454     vSt.RotateZ(ltr->GetAlpha());
2455     vDir.RotateZ(ltr->GetAlpha());
2456
2457     Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2458
2459     if (dist<mindist){
2460       mindist=dist;
2461       itrackMin=itrack;
2462     }
2463
2464   }
2465   itrackMin2=-1;
2466   Float_t mindist2=10;
2467   for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2468     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);  //get the track
2469     if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2470
2471     Double_t deltaZ=ltr->GetZ()-peakpos[2];
2472     if (TMath::Abs(deltaZ)>40) continue;
2473
2474     Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2475     if (dist>1) continue;
2476
2477     if (dist<mindist2){
2478       mindist2=dist;
2479       itrackMin2=itrack;
2480     }
2481   }
2482   mindist=mindist2;
2483   return itrackMin2;
2484
2485 }
2486
2487 //_____________________________________________________________________
2488 Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2489 {
2490   //
2491   // return true if pad is on the edge of a row
2492   //
2493   Int_t edge1   = 0;
2494   if ( pad == edge1 ) return kTRUE;
2495   Int_t edge2   = fROC->GetNPads(sector,row)-1;
2496   if ( pad == edge2 ) return kTRUE;
2497
2498   return kFALSE;
2499 }
2500
2501 //_____________________________________________________________________
2502 TObjArray* AliTPCCalibCE::SetupMeasured()
2503 {
2504   //
2505   // setup array of measured laser tracks and CE
2506   //
2507
2508   TObjArray *arrIdeal    = AliTPCLaserTrack::GetTracks();
2509   TObjArray *arrMeasured = new TObjArray(338);
2510   arrMeasured->SetOwner();
2511   for(Int_t itrack=0;itrack<336;itrack++){
2512     arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2513   }
2514
2515   //"tracks" for CE
2516   AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2517   ltrce->SetId(336);
2518   ltrce->SetSide(0);
2519   ltrce->fVecSec=new TVectorD(557568/2);
2520   ltrce->fVecP2=new TVectorD(557568/2);
2521   ltrce->fVecPhi=new TVectorD(557568/2);
2522   ltrce->fVecGX=new TVectorD(557568/2);
2523   ltrce->fVecGY=new TVectorD(557568/2);
2524   ltrce->fVecGZ=new TVectorD(557568/2);
2525   ltrce->fVecLX=new TVectorD(557568/2);
2526   ltrce->fVecLY=new TVectorD(557568/2);
2527   ltrce->fVecLZ=new TVectorD(557568/2);
2528
2529   arrMeasured->AddAt(ltrce,336); //CE A-Side
2530
2531   ltrce=new AliTPCLaserTrack;
2532   ltrce->SetId(337);
2533   ltrce->SetSide(1);
2534   ltrce->fVecSec=new TVectorD(557568/2);
2535   ltrce->fVecP2=new TVectorD(557568/2);
2536   ltrce->fVecPhi=new TVectorD(557568/2);
2537   ltrce->fVecGX=new TVectorD(557568/2);
2538   ltrce->fVecGY=new TVectorD(557568/2);
2539   ltrce->fVecGZ=new TVectorD(557568/2);
2540   ltrce->fVecLX=new TVectorD(557568/2);
2541   ltrce->fVecLY=new TVectorD(557568/2);
2542   ltrce->fVecLZ=new TVectorD(557568/2);
2543   arrMeasured->AddAt(ltrce,337); //CE C-Side
2544
2545   return arrMeasured;
2546 }
2547
2548 //_____________________________________________________________________
2549 void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2550 {
2551   //
2552   // reset array of measured laser tracks and CE
2553   //
2554   for(Int_t itrack=0;itrack<338;itrack++){
2555     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2556     ltr->fVecSec->Zero();
2557     ltr->fVecP2->Zero();
2558     ltr->fVecPhi->Zero();
2559     ltr->fVecGX->Zero();
2560     ltr->fVecGY->Zero();
2561     ltr->fVecGZ->Zero();
2562     ltr->fVecLX->Zero();
2563     ltr->fVecLY->Zero();
2564     ltr->fVecLZ->Zero();
2565   }
2566 }
2567
2568 //_____________________________________________________________________
2569 void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2570 {
2571   //
2572   // Add ideal CE positions to the ideal track data
2573   //
2574
2575   arr->Expand(338);
2576   //"tracks" for CE
2577   AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2578   ltrceA->SetId(336);
2579   ltrceA->SetSide(0);
2580   ltrceA->fVecSec=new TVectorD(557568/2);
2581   ltrceA->fVecP2=new TVectorD(557568/2);
2582   ltrceA->fVecPhi=new TVectorD(557568/2);
2583   ltrceA->fVecGX=new TVectorD(557568/2);
2584   ltrceA->fVecGY=new TVectorD(557568/2);
2585   ltrceA->fVecGZ=new TVectorD(557568/2);
2586   ltrceA->fVecLX=new TVectorD(557568/2);
2587   ltrceA->fVecLY=new TVectorD(557568/2);
2588   ltrceA->fVecLZ=new TVectorD(557568/2);
2589   arr->AddAt(ltrceA,336); //CE A-Side
2590
2591   AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2592   ltrceC->SetId(337);
2593   ltrceC->SetSide(1);
2594   ltrceC->fVecSec=new TVectorD(557568/2);
2595   ltrceC->fVecP2=new TVectorD(557568/2);
2596   ltrceC->fVecPhi=new TVectorD(557568/2);
2597   ltrceC->fVecGX=new TVectorD(557568/2);
2598   ltrceC->fVecGY=new TVectorD(557568/2);
2599   ltrceC->fVecGZ=new TVectorD(557568/2);
2600   ltrceC->fVecLX=new TVectorD(557568/2);
2601   ltrceC->fVecLY=new TVectorD(557568/2);
2602   ltrceC->fVecLZ=new TVectorD(557568/2);
2603   arr->AddAt(ltrceC,337); //CE C-Side
2604
2605   //Calculate ideal positoins
2606   Float_t gxyz[3];
2607   Float_t lxyz[3];
2608   Int_t channelSideA=0;
2609   Int_t channelSideC=0;
2610   Int_t channelSide=0;
2611   AliTPCLaserTrack *ltrce=0x0;
2612
2613   for (Int_t isector=0; isector<72; ++isector){
2614     Int_t side=((isector/18)%2);
2615     for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2616       for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2617         fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2618         fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2619         if (side==0) {
2620           ltrce=ltrceA;
2621           channelSide=channelSideA;
2622         } else {
2623           ltrce=ltrceC;
2624           channelSide=channelSideC;
2625         }
2626
2627         ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2628         ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2629         ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2630         ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2631         ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2632 //         ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2633         ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2634         ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2635 //         ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2636
2637         if (side==0){
2638           ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2639           ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2640           ++channelSideA;
2641         }
2642         else {
2643           ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2644           ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2645           ++channelSideC;
2646         }
2647       }
2648     }
2649   }
2650
2651
2652 }
2653
2654 //_____________________________________________________________________
2655 void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2656 {
2657   //
2658   // calculate the drift velocity from the reconstructed clusters associated
2659   // to the ideal laser tracks
2660   // use two different fit scenarios: Separate fit for A- and C-Side
2661   //                                  Common fit for A- and C-Side
2662   //
2663
2664   if (!fArrFitGraphs){
2665     fArrFitGraphs=new TObjArray;
2666   }
2667
2668 //   static TLinearFitter fdriftA(5,"hyp4");
2669 //   static TLinearFitter fdriftC(5,"hyp4");
2670 //   static TLinearFitter fdriftAC(6,"hyp5");
2671   Double_t timestamp=fTimeBursts[burst];
2672
2673   static TLinearFitter fdriftA(4,"hyp3");
2674   static TLinearFitter fdriftC(4,"hyp3");
2675   static TLinearFitter fdriftAC(5,"hyp4");
2676   TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2677
2678   Float_t chi2A = 10;
2679   Float_t chi2C = 10;
2680   Float_t chi2AC = 10;
2681   Int_t npointsA=0;
2682   Int_t npointsC=0;
2683   Int_t npointsAC=0;
2684
2685   Double_t minres[3]={20.,1,0.8};
2686   //----
2687   for(Int_t i=0;i<3;i++){
2688
2689     fdriftA.ClearPoints();
2690     fdriftC.ClearPoints();
2691     fdriftAC.ClearPoints();
2692
2693     chi2A = 10;
2694     chi2C = 10;
2695     chi2AC = 10;
2696     npointsA=0;
2697     npointsC=0;
2698     npointsAC=0;
2699
2700     for (Int_t itrack=0; itrack<338; ++itrack){
2701       AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2702       AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2703
2704       //-- Exclude the tracks which has the biggest inclanation angle
2705       if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2706       Int_t clustercounter=0;
2707       Int_t indexMax=159;
2708
2709       //-- exclude the low intensity tracks
2710
2711       for (Int_t index=0; index<indexMax; ++index){
2712
2713         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2714         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2715         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2716
2717         if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2718       }
2719       if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2720       clustercounter=0;
2721
2722
2723       //-- drift length
2724       Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2725
2726       if (itrack>335) indexMax=557568/2;
2727       for (Int_t index=0; index<indexMax; ++index){
2728         Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2729         Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2730         Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2731         Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2732
2733         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2734         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2735         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2736         Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2737
2738       //cut if no track info available
2739         if (iltr->GetBundle()==0) continue;
2740         if (iR<133||mR<133) continue;
2741         if(mltr->fVecP2->GetMatrixArray()[index]>minres[i]) continue;
2742
2743         Double_t ldrift  = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2744         Double_t mdrift  = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2745
2746         //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2747         Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2748
2749         if (iltr->GetSide()==0){
2750           fdriftA.AddPoint(xxx,mdrift,1);
2751         }else{
2752           fdriftC.AddPoint(xxx,mdrift,1);
2753         }
2754 //         Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2755         Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2756         fdriftAC.AddPoint(xxx2,mdrift,1);
2757
2758       }//end index loop
2759     }//end laser track loop
2760
2761   //perform fit
2762     fdriftA.Eval();
2763     fdriftC.Eval();
2764     fdriftAC.Eval();
2765
2766
2767
2768   //get fit values
2769     fdriftA.GetParameters(fitA);
2770     fdriftC.GetParameters(fitC);
2771     fdriftAC.GetParameters(fitAC);
2772
2773   //Parameters:  0 linear offset
2774   //             1 mean drift velocity correction factor
2775   //             2 relative global y gradient
2776   //             3 radial deformation
2777   //             4 IROC/OROC offset
2778
2779 //      FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2780
2781     for (Int_t itrack=0; itrack<338; ++itrack){
2782       AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2783       AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2784
2785       //-- Exclude the tracks which has the biggest inclanation angle
2786       if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2787       Int_t clustercounter=0;
2788       Int_t indexMax=159;
2789
2790       //-- exclude the low intensity tracks
2791
2792       for (Int_t index=0; index<indexMax; ++index){
2793         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2794         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2795         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2796         if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2797       }
2798       if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2799       clustercounter=0;
2800
2801       //-- drift length
2802       Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2803
2804       if (itrack>335) indexMax=557568/2;
2805       for (Int_t index=0; index<indexMax; ++index){
2806         Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2807         Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2808         Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2809         Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2810
2811         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2812         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2813         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2814         Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2815
2816       //cut if no track info available
2817         if (iR<60||mR<60) continue;
2818
2819         Double_t ldrift  = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2820         Double_t mdrift  = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2821
2822         TVectorD *v=&fitA;
2823         if (iltr->GetSide()==1) v=&fitC;
2824 //         Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR)+(*v)[4]*( iltr->fVecSec->GetMatrixArray()[index]>35);
2825         Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2826
2827         mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2828
2829       }
2830     }
2831
2832     fitA.ResizeTo(7);
2833     fitC.ResizeTo(7);
2834     fitAC.ResizeTo(8);
2835
2836 //set statistics values
2837
2838     npointsA= fdriftA.GetNpoints();
2839     if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2840     fitA[5]=npointsA;
2841     fitA[6]=chi2A;
2842
2843     npointsC= fdriftC.GetNpoints();
2844     if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2845     fitC[5]=npointsC;
2846     fitC[6]=chi2C;
2847
2848     npointsAC= fdriftAC.GetNpoints();
2849     if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2850     fitAC[5]=npointsAC;
2851     fitAC[6]=chi2AC;
2852
2853     if (fStreamLevel>2){
2854     //laser track information
2855       (*GetDebugStreamer()) << "DriftV" <<
2856         "iter="   << i <<
2857         "run="    << fRunNumber <<
2858         "time="   << timestamp <<
2859         "fitA.="  << &fitA <<
2860         "fitC.="  << &fitC <<
2861         "fitAC.=" << &fitAC <<
2862         "\n";
2863
2864
2865     }
2866
2867   }
2868 //-----
2869
2870
2871   //Parameters:  0 linear offset (global)
2872   //             1 mean drift velocity correction factor
2873   //             2 relative global y gradient
2874   //             3 radial deformation
2875   //             4 IROC/OROC offset
2876   //             5 linear offset relative A-C
2877
2878   //get graphs
2879   TGraphErrors *grA[7];
2880   TGraphErrors *grC[7];
2881   TGraphErrors *grAC[8];
2882   TString names("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_");
2883   TString namesAC("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_;GRAPH_MEAN_DELAYC_LASER_ALL_");
2884
2885   TObjArray *arrNames=names.Tokenize(";");
2886   TObjArray *arrNamesAC=namesAC.Tokenize(";");
2887
2888   //A-Side graphs
2889   for (Int_t i=0; i<7; ++i){
2890     TString grName=arrNames->UncheckedAt(i)->GetName();
2891     grName+="A";
2892     grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2893     if (!grA[i]){
2894       grA[i]=new TGraphErrors;
2895       grA[i]->SetName(grName.Data());
2896       grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2897       fArrFitGraphs->Add(grA[i]);
2898     }
2899 //     Int_t ipoint=grA[i]->GetN();
2900     Int_t ipoint=burst;
2901     grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2902     grA[i]->SetPointError(ipoint,60,.0001);
2903     if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2904   }
2905
2906   //C-Side graphs
2907   for (Int_t i=0; i<7; ++i){
2908     TString grName=arrNames->UncheckedAt(i)->GetName();
2909     grName+="C";
2910     grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2911     if (!grC[i]){
2912       grC[i]=new TGraphErrors;
2913       grC[i]->SetName(grName.Data());
2914       grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2915       fArrFitGraphs->Add(grC[i]);
2916     }
2917 //     Int_t ipoint=grC[i]->GetN();
2918     Int_t ipoint=burst;
2919     grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2920     grC[i]->SetPointError(ipoint,60,.0001);
2921     if (i<4) grC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2922   }
2923
2924   //AC-Side graphs
2925   for (Int_t i=0; i<8; ++i){
2926     TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2927     grName+="AC";
2928     grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2929     if (!grAC[i]){
2930       grAC[i]=new TGraphErrors;
2931       grAC[i]->SetName(grName.Data());
2932       grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2933       fArrFitGraphs->Add(grAC[i]);
2934     }
2935 //     Int_t ipoint=grAC[i]->GetN();
2936     Int_t ipoint=burst;
2937     grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2938     grAC[i]->SetPointError(ipoint,60,.0001);
2939     if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2940   }
2941
2942   if (fDebugLevel>10){
2943     printf("A side fit parameters:\n");
2944     fitA.Print();
2945     printf("\nC side fit parameters:\n");
2946     fitC.Print();
2947     printf("\nAC side fit parameters:\n");
2948     fitAC.Print();
2949   }
2950   delete arrNames;
2951   delete arrNamesAC;
2952 }
2953
2954 //_____________________________________________________________________
2955 Double_t AliTPCCalibCE::SetBurstHnDrift()
2956 {
2957   //
2958   // Create a new THnSparse for the current burst
2959   // return the time of the current burst
2960   //
2961   Int_t i=0;
2962   for(i=0; i<fTimeBursts.GetNrows(); ++i){
2963     if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2964     if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2965       fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2966       return fTimeBursts(i);
2967     }
2968   }
2969
2970   CreateDVhist();
2971   fArrHnDrift.AddAt(fHnDrift,i);
2972   fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
2973 //   fPeaks[4]=0;
2974   return fTimeStamp;
2975 }
2976
2977 //_____________________________________________________________________
2978 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
2979 {
2980   //
2981   //  Write class to file
2982   //  option can be specified in the dir option:
2983   //  options:
2984   //    name=<objname>: the name of the calibration object in file will be <objname>
2985   //    type=<type>:    the saving type:
2986   //                    0 - write the complte object
2987   //                    1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
2988   //                    2 - like 2, but in addition delete objects that will most probably not be used for calibration
2989   //                    3 - store only calibration output, don't store the reference histograms
2990   //                        and THnSparse (requires Analyse called before)
2991   //
2992   //  NOTE: to read the object back, the ReadFromFile function should be used
2993   //
2994
2995   TString sDir(dir);
2996   TString objName=GetName();
2997   Int_t type=0;
2998
2999   //get options
3000   TObjArray *arr=sDir.Tokenize(",");
3001   TIter next(arr);
3002   TObjString *s;
3003   while ( (s=(TObjString*)next()) ){
3004     TString optString=s->GetString();
3005     optString.Remove(TString::kBoth,' ');
3006     if (optString.BeginsWith("name=")){
3007       objName=optString.ReplaceAll("name=","");
3008     }
3009     if (optString.BeginsWith("type=")){
3010       optString.ReplaceAll("type=","");
3011       type=optString.Atoi();
3012     }
3013   }
3014
3015   if (type==1||type==2) {
3016     //delete most probably not needed stuff
3017     //This requires Analyse to be called after reading the object from file
3018     fCalRocArrayT0.Delete();
3019     fCalRocArrayT0Err.Delete();
3020     fCalRocArrayQ.Delete();
3021     fCalRocArrayRMS.Delete();
3022     fCalRocArrayOutliers.Delete();
3023   }
3024   if (type==2||type==3){
3025     fParamArrayEventPol1.Delete();
3026     fParamArrayEventPol2.Delete();
3027   }
3028
3029   TObjArray histoQArray(72);
3030   TObjArray histoT0Array(72);
3031   TObjArray histoRMSArray(72);
3032   TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3033
3034   //save all large 2D histograms in separte pointers
3035   //to have a smaller memory print when saving the object
3036   if (type==1||type==2||type==3){
3037     for (Int_t i=0; i<72; ++i){
3038       histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3039       histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3040       histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3041     }
3042     fHistoQArray.SetOwner(kFALSE);
3043     fHistoT0Array.SetOwner(kFALSE);
3044     fHistoRMSArray.SetOwner(kFALSE);
3045     fHistoQArray.Clear();
3046     fHistoT0Array.Clear();
3047     fHistoRMSArray.Clear();
3048
3049     for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3050       arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3051     }
3052     fArrHnDrift.SetOwner(kFALSE);
3053     fArrHnDrift.Clear();
3054   }
3055
3056
3057   TDirectory *backup = gDirectory;
3058
3059   TFile f(filename,"recreate");
3060   this->Write(objName.Data());
3061   if (type==1||type==2) {
3062     histoQArray.Write("histoQArray",TObject::kSingleKey);
3063     histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3064     histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3065     arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3066   }
3067
3068   f.Save();
3069   f.Close();
3070
3071   //move histograms back to the object
3072   if (type==1||type==2){
3073     for (Int_t i=0; i<72; ++i){
3074       fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3075       fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3076       fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3077     }
3078     fHistoQArray.SetOwner(kTRUE);
3079     fHistoT0Array.SetOwner(kTRUE);
3080     fHistoRMSArray.SetOwner(kTRUE);
3081
3082     for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3083       fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3084     }
3085     fArrHnDrift.SetOwner(kTRUE);
3086   }
3087
3088   if ( backup ) backup->cd();
3089 }
3090 //_____________________________________________________________________
3091 AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3092 {
3093   //
3094   // Read object from file
3095   // Handle properly if the histogram arrays were stored separately
3096   // call Analyse to make sure to have the calibration relevant information in the object
3097   //
3098
3099   TFile f(filename);
3100   if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3101   TList *l=f.GetListOfKeys();
3102   TIter next(l);
3103   TKey *key=0x0;
3104   TObject *o=0x0;
3105
3106   AliTPCCalibCE *ce=0x0;
3107   TObjArray *histoQArray=0x0;
3108   TObjArray *histoT0Array=0x0;
3109   TObjArray *histoRMSArray=0x0;
3110   TObjArray *arrHnDrift=0x0;
3111
3112   while ( (key=(TKey*)next()) ){
3113     o=key->ReadObj();
3114     if ( o->IsA()==AliTPCCalibCE::Class() ){
3115       ce=(AliTPCCalibCE*)o;
3116     } else if ( o->IsA()==TObjArray::Class() ){
3117       TString name=key->GetName();
3118       if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3119       if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3120       if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3121       if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3122     }
3123   }
3124
3125   if (ce){
3126   //move histograms back to the object
3127     TH1* hist=0x0;
3128     if (histoQArray){
3129       for (Int_t i=0; i<72; ++i){
3130         hist=(TH1*)histoQArray->UncheckedAt(i);
3131         if (hist) hist->SetDirectory(0x0);
3132         ce->fHistoQArray.AddAt(hist,i);
3133       }
3134       ce->fHistoQArray.SetOwner(kTRUE);
3135     }
3136
3137     if (histoT0Array) {
3138       for (Int_t i=0; i<72; ++i){
3139         hist=(TH1*)histoT0Array->UncheckedAt(i);
3140         if (hist) hist->SetDirectory(0x0);
3141         ce->fHistoT0Array.AddAt(hist,i);
3142       }
3143       ce->fHistoT0Array.SetOwner(kTRUE);
3144     }
3145
3146     if (histoRMSArray){
3147       for (Int_t i=0; i<72; ++i){
3148         hist=(TH1*)histoRMSArray->UncheckedAt(i);
3149         if (hist) hist->SetDirectory(0x0);
3150         ce->fHistoRMSArray.AddAt(hist,i);
3151       }
3152       ce->fHistoRMSArray.SetOwner(kTRUE);
3153     }
3154
3155     if (arrHnDrift){
3156       for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3157         THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3158         ce->fArrHnDrift.AddAt(hSparse,i);
3159       }
3160     }
3161
3162     ce->Analyse();
3163   }
3164   f.Close();
3165
3166   return ce;
3167 }