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     Char_t name[255], title[255];
1383
1384     sprintf(name,"hCalib%s%.2d",type,sector);
1385     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1386
1387     // new histogram with Q calib information. One value for each pad!
1388     TH2S* hist = new TH2S(name,title,
1389                           nbinsY, ymin, ymax,
1390                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1391     hist->SetDirectory(0);
1392     arr->AddAt(hist,sector);
1393     return hist;
1394 }
1395 //_____________________________________________________________________
1396 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1397 {
1398     //
1399     // return pointer to T0 histogram
1400     // if force is true create a new histogram if it doesn't exist allready
1401     //
1402     TObjArray *arr = &fHistoT0Array;
1403     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1404 }
1405 //_____________________________________________________________________
1406 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1407 {
1408     //
1409     // return pointer to Q histogram
1410     // if force is true create a new histogram if it doesn't exist allready
1411     //
1412     TObjArray *arr = &fHistoQArray;
1413     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1414 }
1415 //_____________________________________________________________________
1416 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1417 {
1418     //
1419     // return pointer to Q histogram
1420     // if force is true create a new histogram if it doesn't exist allready
1421     //
1422     TObjArray *arr = &fHistoRMSArray;
1423     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1424 }
1425 //_____________________________________________________________________
1426 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1427                               const Char_t *type, Bool_t force)
1428 {
1429     //
1430     // return pointer to TH1S histogram
1431     // if force is true create a new histogram if it doesn't exist allready
1432     //
1433     if ( !force || arr->UncheckedAt(sector) )
1434         return (TH1S*)arr->UncheckedAt(sector);
1435
1436     // if we are forced and histogram doesn't yes exist create it
1437     Char_t name[255], title[255];
1438
1439     sprintf(name,"hCalib%s%.2d",type,sector);
1440     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1441
1442     // new histogram with calib information. One value for each pad!
1443     TH1S* hist = new TH1S(name,title,
1444                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1445     hist->SetDirectory(0);
1446     arr->AddAt(hist,sector);
1447     return hist;
1448 }
1449 //_____________________________________________________________________
1450 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1451 {
1452     //
1453     // return pointer to Q histogram
1454     // if force is true create a new histogram if it doesn't exist allready
1455     //
1456     TObjArray *arr = &fHistoTmean;
1457     return GetHisto(sector, arr, "LastTmean", force);
1458 }
1459 //_____________________________________________________________________
1460 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1461 {
1462     //
1463     // return pointer to Pad Info from 'arr' for the current event and sector
1464     // if force is true create it if it doesn't exist allready
1465     //
1466     if ( !force || arr->UncheckedAt(sector) )
1467         return (TVectorF*)arr->UncheckedAt(sector);
1468
1469     TVectorF *vect = new TVectorF(size);
1470     arr->AddAt(vect,sector);
1471     return vect;
1472 }
1473 //_____________________________________________________________________
1474 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1475 {
1476     //
1477     // return pointer to Pad Times Array for the current event and sector
1478     // if force is true create it if it doesn't exist allready
1479     //
1480     TObjArray *arr = &fPadTimesArrayEvent;
1481     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1482 }
1483 //_____________________________________________________________________
1484 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1485 {
1486     //
1487     // return pointer to Pad Q Array for the current event and sector
1488     // if force is true create it if it doesn't exist allready
1489     // for debugging purposes only
1490     //
1491
1492     TObjArray *arr = &fPadQArrayEvent;
1493     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1494 }
1495 //_____________________________________________________________________
1496 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1497 {
1498     //
1499     // return pointer to Pad RMS Array for the current event and sector
1500     // if force is true create it if it doesn't exist allready
1501     // for debugging purposes only
1502     //
1503     TObjArray *arr = &fPadRMSArrayEvent;
1504     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1505 }
1506 //_____________________________________________________________________
1507 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1508 {
1509     //
1510     // return pointer to Pad RMS Array for the current event and sector
1511     // if force is true create it if it doesn't exist allready
1512     // for debugging purposes only
1513     //
1514     TObjArray *arr = &fPadPedestalArrayEvent;
1515     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1516 }
1517 //_____________________________________________________________________
1518 TVectorF* AliTPCCalibCE::GetTMeanEvents(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 = &fTMeanArrayEvent;
1525     return GetVectSector(sector,arr,100,force);
1526 }
1527 //_____________________________________________________________________
1528 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1529 {
1530     //
1531     // return pointer to the EbyE info of the mean arrival time for 'sector'
1532     // if force is true create it if it doesn't exist allready
1533     //
1534     TObjArray *arr = &fQMeanArrayEvent;
1535     return GetVectSector(sector,arr,100,force);
1536 }
1537 //_____________________________________________________________________
1538 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1539 {
1540     //
1541     // return pointer to ROC Calibration
1542     // if force is true create a new histogram if it doesn't exist allready
1543     //
1544     if ( !force || arr->UncheckedAt(sector) )
1545         return (AliTPCCalROC*)arr->UncheckedAt(sector);
1546
1547     // if we are forced and histogram doesn't yes exist create it
1548
1549     // new AliTPCCalROC for T0 information. One value for each pad!
1550     AliTPCCalROC *croc = new AliTPCCalROC(sector);
1551     arr->AddAt(croc,sector);
1552     return croc;
1553 }
1554 //_____________________________________________________________________
1555 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1556 {
1557     //
1558     // return pointer to Time 0 ROC Calibration
1559     // if force is true create a new histogram if it doesn't exist allready
1560     //
1561     TObjArray *arr = &fCalRocArrayT0;
1562     return GetCalRoc(sector, arr, force);
1563 }
1564 //_____________________________________________________________________
1565 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1566 {
1567     //
1568     // return pointer to the error of Time 0 ROC Calibration
1569     // if force is true create a new histogram if it doesn't exist allready
1570     //
1571     TObjArray *arr = &fCalRocArrayT0Err;
1572     return GetCalRoc(sector, arr, force);
1573 }
1574 //_____________________________________________________________________
1575 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1576 {
1577     //
1578     // return pointer to T0 ROC Calibration
1579     // if force is true create a new histogram if it doesn't exist allready
1580     //
1581     TObjArray *arr = &fCalRocArrayQ;
1582     return GetCalRoc(sector, arr, force);
1583 }
1584 //_____________________________________________________________________
1585 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1586 {
1587     //
1588     // return pointer to signal width ROC Calibration
1589     // if force is true create a new histogram if it doesn't exist allready
1590     //
1591     TObjArray *arr = &fCalRocArrayRMS;
1592     return GetCalRoc(sector, arr, force);
1593 }
1594 //_____________________________________________________________________
1595 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1596 {
1597     //
1598     // return pointer to Outliers
1599     // if force is true create a new histogram if it doesn't exist allready
1600     //
1601     TObjArray *arr = &fCalRocArrayOutliers;
1602     return GetCalRoc(sector, arr, force);
1603 }
1604 //_____________________________________________________________________
1605 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1606 {
1607     //
1608     // return pointer to TObjArray of fit parameters
1609     // if force is true create a new histogram if it doesn't exist allready
1610     //
1611     if ( !force || arr->UncheckedAt(sector) )
1612         return (TObjArray*)arr->UncheckedAt(sector);
1613
1614     // if we are forced and array doesn't yes exist create it
1615
1616     // new TObjArray for parameters
1617     TObjArray *newArr = new TObjArray;
1618     arr->AddAt(newArr,sector);
1619     return newArr;
1620 }
1621 //_____________________________________________________________________
1622 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1623 {
1624     //
1625     // return pointer to TObjArray of fit parameters from plane fit
1626     // if force is true create a new histogram if it doesn't exist allready
1627     //
1628     TObjArray *arr = &fParamArrayEventPol1;
1629     return GetParamArray(sector, arr, force);
1630 }
1631 //_____________________________________________________________________
1632 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1633 {
1634     //
1635     // return pointer to TObjArray of fit parameters from parabola fit
1636     // if force is true create a new histogram if it doesn't exist allready
1637     //
1638     TObjArray *arr = &fParamArrayEventPol2;
1639     return GetParamArray(sector, arr, force);
1640 }
1641
1642 //_____________________________________________________________________
1643 void AliTPCCalibCE::CreateDVhist()
1644 {
1645   //
1646   // Setup the THnSparse for the drift velocity determination
1647   //
1648
1649   //HnSparse bins
1650   //roc, row, pad, timebin, timestamp
1651   TTimeStamp begin(2010,01,01,0,0,0);
1652   TTimeStamp end(2035,01,01,0,0,0);
1653   Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1654
1655   Int_t    bins[kHnBinsDV] = { 72,  96,  140,  1030, nbinsTime};
1656   Double_t xmin[kHnBinsDV] = { 0.,  0.,   0.,    0., (Double_t)begin.GetSec()};
1657   Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1658
1659   fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1660   fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1661   fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1662   fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1663   fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1664   fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
1665
1666 }
1667
1668 //_____________________________________________________________________
1669 void AliTPCCalibCE::ResetEvent()
1670 {
1671     //
1672     //  Reset global counters  -- Should be called before each event is processed
1673     //
1674     fLastSector=-1;
1675     fCurrentSector=-1;
1676     fCurrentRow=-1;
1677     fCurrentChannel=-1;
1678
1679     ResetPad();
1680
1681     fPadTimesArrayEvent.Delete();
1682     fPadQArrayEvent.Delete();
1683     fPadRMSArrayEvent.Delete();
1684     fPadPedestalArrayEvent.Delete();
1685
1686     for ( Int_t i=0; i<72; ++i ){
1687         fVTime0Offset.GetMatrixArray()[i]=0;
1688         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1689         fVMeanQ.GetMatrixArray()[i]=0;
1690         fVMeanQCounter.GetMatrixArray()[i]=0;
1691     }
1692 }
1693 //_____________________________________________________________________
1694 void AliTPCCalibCE::ResetPad()
1695 {
1696     //
1697     //  Reset pad infos -- Should be called after a pad has been processed
1698     //
1699     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1700       fPadSignal[i] = 0;
1701     fMaxTimeBin   = -1;
1702     fMaxPadSignal = -1;
1703     fPadPedestal  = -1;
1704     fPadNoise     = -1;
1705 }
1706 //_____________________________________________________________________
1707 void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
1708 {
1709   //
1710   //  Merge ce to the current AliTPCCalibCE
1711   //
1712   MergeBase(ce);
1713   Int_t nCEevents = ce->GetNeventsProcessed();
1714
1715   if (fProcessOld&&ce->fProcessOld){
1716   //merge histograms
1717     for (Int_t iSec=0; iSec<72; ++iSec){
1718       TH2S *hRefQmerge   = ce->GetHistoQ(iSec);
1719       TH2S *hRefT0merge  = ce->GetHistoT0(iSec);
1720       TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1721
1722
1723       if ( hRefQmerge ){
1724         TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1725         TH2S *hRefQ   = GetHistoQ(iSec);
1726         if ( hRefQ ) hRefQ->Add(hRefQmerge);
1727         else {
1728           TH2S *hist = new TH2S(*hRefQmerge);
1729           hist->SetDirectory(0);
1730           fHistoQArray.AddAt(hist, iSec);
1731         }
1732         hRefQmerge->SetDirectory(dir);
1733       }
1734       if ( hRefT0merge ){
1735         TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1736         TH2S *hRefT0  = GetHistoT0(iSec);
1737         if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1738         else {
1739           TH2S *hist = new TH2S(*hRefT0merge);
1740           hist->SetDirectory(0);
1741           fHistoT0Array.AddAt(hist, iSec);
1742         }
1743         hRefT0merge->SetDirectory(dir);
1744       }
1745       if ( hRefRMSmerge ){
1746         TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1747         TH2S *hRefRMS = GetHistoRMS(iSec);
1748         if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1749         else {
1750           TH2S *hist = new TH2S(*hRefRMSmerge);
1751           hist->SetDirectory(0);
1752           fHistoRMSArray.AddAt(hist, iSec);
1753         }
1754         hRefRMSmerge->SetDirectory(dir);
1755       }
1756
1757     }
1758
1759     // merge time information
1760
1761
1762     for (Int_t iSec=0; iSec<72; ++iSec){
1763       TObjArray *arrPol1CE  = ce->GetParamArrayPol1(iSec);
1764       TObjArray *arrPol2CE  = ce->GetParamArrayPol2(iSec);
1765       TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1766       TVectorF *vMeanQCE    = ce->GetQMeanEvents(iSec);
1767
1768       TObjArray *arrPol1  = 0x0;
1769       TObjArray *arrPol2  = 0x0;
1770       TVectorF *vMeanTime = 0x0;
1771       TVectorF *vMeanQ    = 0x0;
1772
1773   //resize arrays
1774       if ( arrPol1CE && arrPol2CE ){
1775         arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1776         arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1777         arrPol1->Expand(fNevents+nCEevents);
1778         arrPol2->Expand(fNevents+nCEevents);
1779       }
1780       if ( vMeanTimeCE && vMeanQCE ){
1781         vMeanTime = GetTMeanEvents(iSec,kTRUE);
1782         vMeanQ    = GetQMeanEvents(iSec,kTRUE);
1783         vMeanTime->ResizeTo(fNevents+nCEevents);
1784         vMeanQ->ResizeTo(fNevents+nCEevents);
1785       }
1786
1787       for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1788         if ( arrPol1CE && arrPol2CE ){
1789           TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1790           TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1791           if ( paramPol1 && paramPol2 ){
1792             GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1793             GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1794           }
1795         }
1796         if ( vMeanTimeCE && vMeanQCE ){
1797           vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1798           vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1799         }
1800       }
1801     }
1802
1803
1804
1805     const TVectorD&  eventTimes  = ce->fVEventTime;
1806     const TVectorD&  eventIds    = ce->fVEventNumber;
1807     const TVectorF&  time0SideA  = ce->fVTime0SideA;
1808     const TVectorF&  time0SideC  = ce->fVTime0SideC;
1809     fVEventTime.ResizeTo(fNevents+nCEevents);
1810     fVEventNumber.ResizeTo(fNevents+nCEevents);
1811     fVTime0SideA.ResizeTo(fNevents+nCEevents);
1812     fVTime0SideC.ResizeTo(fNevents+nCEevents);
1813
1814     for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1815       Double_t evTime     = eventTimes.GetMatrixArray()[iEvent];
1816       Double_t evId       = eventIds.GetMatrixArray()[iEvent];
1817       Float_t  t0SideA    = time0SideA.GetMatrixArray()[iEvent];
1818       Float_t  t0SideC    = time0SideC.GetMatrixArray()[iEvent];
1819
1820       fVEventTime.GetMatrixArray()[fNevents+iEvent]   = evTime;
1821       fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1822       fVTime0SideA.GetMatrixArray()[fNevents+iEvent]  = t0SideA;
1823       fVTime0SideC.GetMatrixArray()[fNevents+iEvent]  = t0SideC;
1824     }
1825   }
1826
1827   if (fProcessNew&&ce->fProcessNew) {
1828     if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1829       AliError("Number of bursts in the instances to merge are different. No merging done!");
1830     } else {
1831       for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1832         THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1833         THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1834         if (h && hce) h->Add(hce);
1835         else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
1836       }
1837       //TODO: What to do with fTimeBursts???
1838     }
1839   }
1840
1841   fNevents+=nCEevents; //increase event counter
1842 }
1843
1844 //_____________________________________________________________________
1845 Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1846 {
1847   //
1848   // Merge all objects of this type in list
1849   //
1850
1851   Long64_t nmerged=1;
1852
1853   TIter next(list);
1854   AliTPCCalibCE *ce=0;
1855   TObject *o=0;
1856
1857   while ( (o=next()) ){
1858     ce=dynamic_cast<AliTPCCalibCE*>(o);
1859     if (ce){
1860       Merge(ce);
1861       ++nmerged;
1862     }
1863   }
1864
1865   return nmerged;
1866 }
1867
1868 //_____________________________________________________________________
1869 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1870 {
1871   //
1872   // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1873   // or side (-1: A-Side, -2: C-Side)
1874   // xVariable:    0-event time, 1-event id, 2-internal event counter
1875   // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1876   // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1877   //                               0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1878   //                               not used for mean time and mean Q )
1879   // for an example see class description at the beginning
1880   //
1881
1882   Double_t *x = new Double_t[fNevents];
1883   Double_t *y = new Double_t[fNevents];
1884
1885   TVectorD *xVar = 0x0;
1886   TObjArray *aType = 0x0;
1887   Int_t npoints=0;
1888
1889   // sanity checks
1890   if ( (sector<-2)   || (sector>71)   ) return 0x0;  //sector outside valid range
1891   if ( (xVariable<0) || (xVariable>2) ) return 0x0;  //invalid x-variable
1892   if ( (fitType<0)   || (fitType>3)   ) return 0x0;  //invalid fit type
1893   if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1894   if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1895   if ( sector<0 && fitType!=2) return 0x0;  //for side wise information only mean time is available
1896
1897   if (sector>=0){
1898     if ( fitType==0 ){
1899       if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1900       aType = &fParamArrayEventPol1;
1901       if ( aType->At(sector)==0x0 ) return 0x0;
1902     }
1903     else if ( fitType==1 ){
1904       if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1905       aType = &fParamArrayEventPol2;
1906       if ( aType->At(sector)==0x0 ) return 0x0;
1907     }
1908
1909   }
1910   if ( xVariable == 0 ) xVar = &fVEventTime;
1911   if ( xVariable == 1 ) xVar = &fVEventNumber;
1912   if ( xVariable == 2 ) {
1913     xVar = new TVectorD(fNevents);
1914     for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1915   }
1916
1917   for (Int_t ievent =0; ievent<fNevents; ++ievent){
1918     if ( fitType<2 ){
1919       TObjArray *events = (TObjArray*)(aType->At(sector));
1920       if ( events->GetSize()<=ievent ) break;
1921       TVectorD *v = (TVectorD*)(events->At(ievent));
1922       if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1923     } else if (fitType == 2) {
1924       Double_t xValue=(*xVar)[ievent];
1925       Double_t yValue=0;
1926       if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1927       else if (sector==-1) yValue=fVTime0SideA(ievent);
1928       else if (sector==-2) yValue=fVTime0SideC(ievent);
1929       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1930     }else if (fitType == 3) {
1931       Double_t xValue=(*xVar)[ievent];
1932       Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1933       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1934     }
1935   }
1936
1937   TGraph *gr = new TGraph(npoints);
1938     //sort xVariable increasing
1939   Int_t    *sortIndex = new Int_t[npoints];
1940   TMath::Sort(npoints,x,sortIndex);
1941   for (Int_t i=0;i<npoints;++i){
1942     gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1943   }
1944
1945
1946   if ( xVariable == 2 ) delete xVar;
1947   delete [] x;
1948   delete [] y;
1949   delete [] sortIndex;
1950   return gr;
1951 }
1952 //_____________________________________________________________________
1953 void AliTPCCalibCE::Analyse()
1954 {
1955   //
1956   //  Calculate calibration constants
1957   //
1958
1959   if (fProcessOld){
1960     TVectorD paramQ(3);
1961     TVectorD paramT0(3);
1962     TVectorD paramRMS(3);
1963     TMatrixD dummy(3,3);
1964
1965     Float_t channelCounter=0;
1966     fMeanT0rms=0;
1967     fMeanQrms=0;
1968     fMeanRMSrms=0;
1969
1970     for (Int_t iSec=0; iSec<72; ++iSec){
1971       TH2S *hT0 = GetHistoT0(iSec);
1972       if (!hT0 ) continue;
1973
1974       AliTPCCalROC *rocQ     = GetCalRocQ  (iSec,kTRUE);
1975       AliTPCCalROC *rocT0    = GetCalRocT0 (iSec,kTRUE);
1976       AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1977       AliTPCCalROC *rocRMS   = GetCalRocRMS(iSec,kTRUE);
1978       AliTPCCalROC *rocOut   = GetCalRocOutliers(iSec,kTRUE);
1979
1980       TH2S *hQ   = GetHistoQ(iSec);
1981       TH2S *hRMS = GetHistoRMS(iSec);
1982
1983       Short_t *arrayhQ   = hQ->GetArray();
1984       Short_t *arrayhT0  = hT0->GetArray();
1985       Short_t *arrayhRMS = hRMS->GetArray();
1986
1987       UInt_t nChannels = fROC->GetNChannels(iSec);
1988
1989   //debug
1990       Int_t row=0;
1991       Int_t pad=0;
1992       Int_t padc=0;
1993   //! debug
1994
1995       for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1996
1997
1998         Float_t cogTime0 = -1000;
1999         Float_t cogQ     = -1000;
2000         Float_t cogRMS   = -1000;
2001         Float_t cogOut   = 0;
2002         Float_t rms      = 0;
2003         Float_t rmsT0    = 0;
2004
2005
2006         Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
2007         Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
2008         Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
2009
2010         cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2011         fMeanQrms+=rms;
2012         cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2013         fMeanT0rms+=rmsT0;
2014         cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2015         fMeanRMSrms+=rms;
2016         channelCounter++;
2017
2018       /*
2019              //outlier specifications
2020          if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2021          cogOut = 1;
2022           cogTime0 = 0;
2023           cogQ     = 0;
2024           cogRMS   = 0;
2025       }
2026       */
2027         rocQ->SetValue(iChannel, cogQ*cogQ);
2028         rocT0->SetValue(iChannel, cogTime0);
2029         rocT0Err->SetValue(iChannel, rmsT0);
2030         rocRMS->SetValue(iChannel, cogRMS);
2031         rocOut->SetValue(iChannel, cogOut);
2032
2033
2034       //debug
2035         if ( GetStreamLevel() > 0 ){
2036           TTreeSRedirector *streamer=GetDebugStreamer();
2037           if ( streamer ) {
2038
2039             while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2040             pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2041             padc = pad-(fROC->GetNPads(iSec,row)/2);
2042
2043             (*streamer) << "DataEnd" <<
2044               "Sector="  << iSec      <<
2045               "Pad="     << pad       <<
2046               "PadC="    << padc      <<
2047               "Row="     << row       <<
2048               "PadSec="  << iChannel   <<
2049               "Q="       << cogQ      <<
2050               "T0="      << cogTime0  <<
2051               "RMS="     << cogRMS    <<
2052               "\n";
2053           }
2054         }
2055       //! debug
2056
2057       }
2058
2059     }
2060     if ( channelCounter>0 ){
2061       fMeanT0rms/=channelCounter;
2062       fMeanQrms/=channelCounter;
2063       fMeanRMSrms/=channelCounter;
2064     }
2065     //   if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2066     //    delete fDebugStreamer;
2067     //    fDebugStreamer = 0x0;
2068     fVEventTime.ResizeTo(fNevents);
2069     fVEventNumber.ResizeTo(fNevents);
2070     fVTime0SideA.ResizeTo(fNevents);
2071     fVTime0SideC.ResizeTo(fNevents);
2072   }
2073
2074   if (fProcessNew&&fAnalyseNew){
2075     AnalyseTrack();
2076     for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2077       THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2078       h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2079     }
2080   }
2081 }
2082
2083
2084
2085
2086 //
2087 // New functions that also use the laser tracks
2088 //
2089
2090
2091
2092 //_____________________________________________________________________
2093 void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2094 {
2095   //
2096   //Find the local maximums for the projections to each axis
2097   //
2098
2099   //find laser layer positoins
2100   fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2101   FindLaserLayers();
2102   THnSparse *hProj=fHnDrift;
2103   Double_t posCE[4]={0.,0.,0.,0.};
2104   Double_t widthCE[4]={0.,0.,0.,0.};
2105
2106 //   if(fPeaks[4]!=0){
2107   // find central electrode position once more, separately for IROC, OROC, A-, C-Side
2108
2109   for (Int_t i=0; i<4; ++i){
2110     hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2111     TH1 *h=fHnDrift->Projection(3);
2112     h->GetXaxis()->SetRangeUser(fPeaks[4]-fPeakWidths[4],fPeaks[4]+fPeakWidths[4]);
2113     Int_t nbinMax=h->GetMaximumBin();
2114     Double_t maximum=h->GetMaximum();
2115 //     Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2116 //     if (nbinMax<700||maximum<maxExpected) continue;
2117     Double_t xbinMax=h->GetBinCenter(nbinMax);
2118     TF1 fgaus("gaus","gaus",xbinMax-10,xbinMax+10);
2119     fgaus.SetParameters(maximum,xbinMax,2);
2120     fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2121     fgaus.SetParLimits(2,0.2,4.);
2122     h->Fit(&fgaus,"RQN");
2123 //     Double_t deltaX=4*fgaus.GetParameter(2);
2124 //     xbinMax=fgaus.GetParameter(1);
2125     delete h;
2126     posCE[i]=fgaus.GetParameter(1);
2127     widthCE[i]=4*fgaus.GetParameter(2);
2128     hProj->GetAxis(0)->SetRangeUser(0,72);
2129   }
2130 //   }
2131   //Current drift velocity
2132   Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2133 //   cout<<"vdrift="<<vdrift<<endl;
2134
2135   AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2136   //loop over all entries in the histogram
2137   Int_t coord[5];
2138   for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2139     //get entry position and content
2140     Double_t adc=hProj->GetBinContent(ichunk,coord);
2141
2142
2143     Int_t sector = coord[0]-1;
2144     Int_t row    = coord[1]-1;
2145     Int_t pad    = coord[2]-1;
2146     Int_t timeBin= coord[3]-1;
2147     Double_t time   = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2148     Int_t side   = (sector/18)%2;
2149 //     return;
2150 //     fPeaks[4]=(UInt_t)posCE[sector/18];
2151 //     fPeakWidths[4]=(UInt_t)widthCE[sector/18];
2152
2153     //cuts
2154     if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2155     if (adc < 5 ) continue;
2156     if (IsEdgePad(sector,row,pad)) continue;
2157 //     if (!IsPeakInRange(timeBin)) continue;
2158 //     if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2159 //     if (isCE&&(sector!=0)) continue;
2160
2161     Int_t padmin=-2, padmax=2;
2162     Int_t timemin=-2, timemax=2;
2163     Int_t minsumperpad=2;
2164     //CE or laser tracks
2165     Bool_t isCE=kFALSE;
2166     if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2167       isCE=kTRUE;
2168       padmin=0;
2169       padmax=0;
2170       timemin=-3;
2171       timemax=7;
2172     }
2173
2174     //
2175     // Find local maximum and cogs
2176     //
2177     Bool_t isMaximum=kTRUE;
2178     Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2179     Double_t cogY=0, rmsY=0;
2180     Int_t npart=0;
2181
2182     // for position calculation use
2183     for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2184       Float_t lxyz[3];
2185       fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2186
2187       for(Int_t itime=timemin;itime<=timemax;++itime){
2188
2189         Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2190         Double_t val=hProj->GetBinContent(a);
2191         totalmass+=val;
2192
2193         tbcm +=(timeBin+itime)*val;
2194         padcm+=(pad+ipad)*val;
2195         cogY +=lxyz[1]*val;
2196
2197         rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2198         rmspad+=(pad+ipad)*(pad+ipad)*val;
2199         rmsY  +=lxyz[1]*lxyz[1]*val;
2200
2201         if (val>0) ++npart;
2202         if (val>adc) {
2203           isMaximum=kFALSE;
2204           break;
2205         }
2206       }
2207       if (!isMaximum)  break;
2208     }
2209
2210     if (!isMaximum||!npart)  continue;
2211     if (totalmass<npart*minsumperpad) continue;
2212     if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2213                                     //for CE we want only one pad by construction
2214
2215     tbcm/=totalmass;
2216     padcm/=totalmass;
2217     cogY/=totalmass;
2218
2219     rmstb/=totalmass;
2220     rmspad/=totalmass;
2221     rmsY/=totalmass;
2222
2223     rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2224     rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2225     rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2226
2227     Int_t cog=TMath::Nint(padcm);
2228
2229     // timebin --> z position
2230     Float_t zlength=fParam->GetZLength(side);
2231 //     Float_t timePos=tbcm+(1000-fPeaks[4])
2232     // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2233     Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2234
2235     // local to global transformation--> x and y positions
2236     Float_t padlxyz[3];
2237     fROC->GetPositionLocal(sector,row,pad,padlxyz);
2238
2239     Double_t gxyz[3]={padlxyz[0],cogY,gz};
2240     Double_t lxyz[3]={padlxyz[0],cogY,gz};
2241     Double_t igxyz[3]={0,0,0};
2242     AliTPCTransform t1;
2243     t1.RotatedGlobal2Global(sector,gxyz);
2244
2245     Double_t mindist=0;
2246     Int_t trackID=-1;
2247     Int_t trackID2=-1;
2248
2249     //find track id and index of the position in the track (row)
2250     Int_t index=0;
2251     if (!isCE){
2252       index=row+(sector>35)*fROC->GetNRows(0);
2253       trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2254     } else {
2255       trackID=336+((sector/18)%2);
2256       index= fROC->GetRowIndexes(sector)[row]+pad; //  global pad position in sector
2257       if (sector<36) {
2258         index+=(sector%18)*fROC->GetNChannels(sector);
2259       } else {
2260         index+=18*fROC->GetNChannels(0);
2261         index+=(sector%18)*fROC->GetNChannels(sector);
2262       }
2263       //TODO: find out about the multiple peaks in the CE
2264 //       mindist=TMath::Abs(fPeaks[4]-tbcm);
2265       mindist=1.;
2266     }
2267
2268     // fill track vectors
2269     if (trackID>0){
2270       AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2271       Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2272
2273
2274 //      travel time effect of light includes
2275
2276       Double_t raylength=ltr->GetRayLength();
2277       Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2278       ltr->GetXYZ(globmir);
2279       if(trackID<336){
2280         if(side==0){
2281           gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2282                                        +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2283                                        +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2284         }
2285         else {
2286           gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2287                                        +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2288                                        +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2289         }
2290       }
2291
2292       if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2293         ltr->fVecSec->GetMatrixArray()[index]=sector;
2294         ltr->fVecP2->GetMatrixArray()[index]=0;
2295         ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2296         ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2297         ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2298         ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2299         ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2300         ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2301         ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2302 //         ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2303       }
2304       TObjArray *arr=AliTPCLaserTrack::GetTracks();
2305       ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2306       igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2307       igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2308       igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2309     }
2310
2311
2312     if (fStreamLevel>4){
2313       (*GetDebugStreamer()) << "clusters" <<
2314         "run="   << fRunNumber <<
2315         "timestamp=" << timestamp <<
2316         "burst="     << burst     <<
2317         "side="      << side      <<
2318         "sec="       << sector    <<
2319         "row="       << row       <<
2320         "pad="       << pad       <<
2321         "padCog="    << cog       <<
2322         "timebin="   << timeBin   <<
2323         "cogCE="     << posCE[sector/18] <<
2324         "withCE="    << widthCE[sector/18] <<
2325         "index="     << index     <<
2326
2327         "padcm="     << padcm     <<
2328         "rmspad="    << rmspad    <<
2329
2330         "cogtb="     << tbcm      <<
2331         "rmstb="     << rmstb     <<
2332
2333         "npad="      << npart     <<
2334
2335         "lx="        << padlxyz[0]<<
2336         "ly="        << cogY      <<
2337         "lypad="     << padlxyz[1]<<
2338         "rmsY="      << rmsY      <<
2339
2340         "gx="        << gxyz[0]   <<
2341         "gy="        << gxyz[1]   <<
2342         "gz="        << gxyz[2]   <<
2343
2344         "igx="        << igxyz[0] <<
2345         "igy="        << igxyz[1] <<
2346         "igz="        << igxyz[2] <<
2347
2348         "mind="      << mindist   <<
2349         "max="       << adc       <<
2350         "trackid="   << trackID   <<
2351         "trackid2="   << trackID2   <<
2352         "npart="     << npart     <<
2353         "\n";
2354     } // end stream levelmgz.fElements
2355
2356   }
2357
2358 }
2359
2360 //_____________________________________________________________________
2361 void AliTPCCalibCE::AnalyseTrack()
2362 {
2363   //
2364   //  Analyse the tracks
2365   //
2366
2367
2368   AliTPCLaserTrack::LoadTracks();
2369 //   AliTPCParam *param=0x0;
2370 //   //cdb run number
2371 //   AliCDBManager *man=AliCDBManager::Instance();
2372 //   if (man->GetDefaultStorage()){
2373 //     AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2374 //     if (entry){
2375 //       entry->SetOwner(kTRUE);
2376 //       param = (AliTPCParam*)(entry->GetObject()->Clone());
2377 //     }
2378 //   }
2379 //   if (param){
2380 //     if (fParam) delete fParam;
2381 //     fParam=param;
2382 //   } else {
2383 //     AliError("Could not get updated AliTPCParam from OCDB!!!");
2384 //   }
2385
2386   //Measured and ideal laser tracks
2387   TObjArray* arrMeasured = SetupMeasured();
2388   TObjArray* arrIdeal    = AliTPCLaserTrack::GetTracks();
2389   AddCEtoIdeal(arrIdeal);
2390
2391   //find bursts and loop over them
2392   for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2393     Double_t timestamp=fTimeBursts[iburst];
2394     AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2395     fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2396     if (!fHnDrift) continue;
2397     UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2398     if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2399     fBinsLastAna[iburst]=entries;
2400
2401     for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2402 //     if (iburst==0) FindLaserLayers();
2403
2404     //reset laser tracks
2405     ResetMeasured(arrMeasured);
2406
2407     //find clusters and associate to the tracks
2408     FindLocalMaxima(arrMeasured, timestamp, iburst);
2409
2410     //calculate drift velocity
2411     CalculateDV(arrIdeal,arrMeasured,iburst);
2412
2413     //Dump information to file if requested
2414     if (fStreamLevel>2){
2415       printf("make tree\n");
2416       //laser track information
2417
2418       for (Int_t itrack=0; itrack<338; ++itrack){
2419         TObject *iltr=arrIdeal->UncheckedAt(itrack);
2420         TObject *mltr=arrMeasured->UncheckedAt(itrack);
2421         (*GetDebugStreamer()) << "tracks" <<
2422           "run="   << fRunNumber <<
2423           "time=" << timestamp <<
2424           "burst="<< iburst <<
2425           "iltr.=" << iltr <<
2426           "mltr.=" << mltr <<
2427           "\n";
2428       }
2429     }
2430   }
2431   if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2432 }
2433
2434 //_____________________________________________________________________
2435 Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2436                                       const Double_t *peakposloc, Int_t &itrackMin2)
2437 {
2438   //
2439   //  Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2440   //
2441
2442
2443   TObjArray *arr=AliTPCLaserTrack::GetTracks();
2444   TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2445   TVector3 vDir;
2446   TVector3 vSt;
2447
2448   Int_t firstbeam=0;
2449   Int_t lastbeam=336/2;
2450   if ( (sector/18)%2 ) {
2451     firstbeam=336/2;
2452     lastbeam=336;
2453   }
2454
2455   mindist=1000000;
2456   Int_t itrackMin=-1;
2457   for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2458     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);  //get the track
2459 //     if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2460     vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2461     Double_t deltaZ=ltr->GetZ()-peakpos[2];
2462     if (TMath::Abs(deltaZ)>40) continue;
2463     vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2464     vSt.RotateZ(ltr->GetAlpha());
2465     vDir.RotateZ(ltr->GetAlpha());
2466
2467     Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2468
2469     if (dist<mindist){
2470       mindist=dist;
2471       itrackMin=itrack;
2472     }
2473
2474   }
2475   itrackMin2=-1;
2476   Float_t mindist2=10;
2477   for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2478     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);  //get the track
2479     if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2480
2481     Double_t deltaZ=ltr->GetZ()-peakpos[2];
2482     if (TMath::Abs(deltaZ)>40) continue;
2483
2484     Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2485     if (dist>1) continue;
2486
2487     if (dist<mindist2){
2488       mindist2=dist;
2489       itrackMin2=itrack;
2490     }
2491   }
2492   mindist=mindist2;
2493   return itrackMin2;
2494
2495 }
2496
2497 //_____________________________________________________________________
2498 Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2499 {
2500   //
2501   // return true if pad is on the edge of a row
2502   //
2503   Int_t edge1   = 0;
2504   if ( pad == edge1 ) return kTRUE;
2505   Int_t edge2   = fROC->GetNPads(sector,row)-1;
2506   if ( pad == edge2 ) return kTRUE;
2507
2508   return kFALSE;
2509 }
2510
2511 //_____________________________________________________________________
2512 TObjArray* AliTPCCalibCE::SetupMeasured()
2513 {
2514   //
2515   // setup array of measured laser tracks and CE
2516   //
2517
2518   TObjArray *arrIdeal    = AliTPCLaserTrack::GetTracks();
2519   TObjArray *arrMeasured = new TObjArray(338);
2520   arrMeasured->SetOwner();
2521   for(Int_t itrack=0;itrack<336;itrack++){
2522     arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2523   }
2524
2525   //"tracks" for CE
2526   AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2527   ltrce->SetId(336);
2528   ltrce->SetSide(0);
2529   ltrce->fVecSec=new TVectorD(557568/2);
2530   ltrce->fVecP2=new TVectorD(557568/2);
2531   ltrce->fVecPhi=new TVectorD(557568/2);
2532   ltrce->fVecGX=new TVectorD(557568/2);
2533   ltrce->fVecGY=new TVectorD(557568/2);
2534   ltrce->fVecGZ=new TVectorD(557568/2);
2535   ltrce->fVecLX=new TVectorD(557568/2);
2536   ltrce->fVecLY=new TVectorD(557568/2);
2537   ltrce->fVecLZ=new TVectorD(557568/2);
2538
2539   arrMeasured->AddAt(ltrce,336); //CE A-Side
2540
2541   ltrce=new AliTPCLaserTrack;
2542   ltrce->SetId(337);
2543   ltrce->SetSide(1);
2544   ltrce->fVecSec=new TVectorD(557568/2);
2545   ltrce->fVecP2=new TVectorD(557568/2);
2546   ltrce->fVecPhi=new TVectorD(557568/2);
2547   ltrce->fVecGX=new TVectorD(557568/2);
2548   ltrce->fVecGY=new TVectorD(557568/2);
2549   ltrce->fVecGZ=new TVectorD(557568/2);
2550   ltrce->fVecLX=new TVectorD(557568/2);
2551   ltrce->fVecLY=new TVectorD(557568/2);
2552   ltrce->fVecLZ=new TVectorD(557568/2);
2553   arrMeasured->AddAt(ltrce,337); //CE C-Side
2554
2555   return arrMeasured;
2556 }
2557
2558 //_____________________________________________________________________
2559 void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2560 {
2561   //
2562   // reset array of measured laser tracks and CE
2563   //
2564   for(Int_t itrack=0;itrack<338;itrack++){
2565     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2566     ltr->fVecSec->Zero();
2567     ltr->fVecP2->Zero();
2568     ltr->fVecPhi->Zero();
2569     ltr->fVecGX->Zero();
2570     ltr->fVecGY->Zero();
2571     ltr->fVecGZ->Zero();
2572     ltr->fVecLX->Zero();
2573     ltr->fVecLY->Zero();
2574     ltr->fVecLZ->Zero();
2575   }
2576 }
2577
2578 //_____________________________________________________________________
2579 void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2580 {
2581   //
2582   // Add ideal CE positions to the ideal track data
2583   //
2584
2585   arr->Expand(338);
2586   //"tracks" for CE
2587   AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2588   ltrceA->SetId(336);
2589   ltrceA->SetSide(0);
2590   ltrceA->fVecSec=new TVectorD(557568/2);
2591   ltrceA->fVecP2=new TVectorD(557568/2);
2592   ltrceA->fVecPhi=new TVectorD(557568/2);
2593   ltrceA->fVecGX=new TVectorD(557568/2);
2594   ltrceA->fVecGY=new TVectorD(557568/2);
2595   ltrceA->fVecGZ=new TVectorD(557568/2);
2596   ltrceA->fVecLX=new TVectorD(557568/2);
2597   ltrceA->fVecLY=new TVectorD(557568/2);
2598   ltrceA->fVecLZ=new TVectorD(557568/2);
2599   arr->AddAt(ltrceA,336); //CE A-Side
2600
2601   AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2602   ltrceC->SetId(337);
2603   ltrceC->SetSide(1);
2604   ltrceC->fVecSec=new TVectorD(557568/2);
2605   ltrceC->fVecP2=new TVectorD(557568/2);
2606   ltrceC->fVecPhi=new TVectorD(557568/2);
2607   ltrceC->fVecGX=new TVectorD(557568/2);
2608   ltrceC->fVecGY=new TVectorD(557568/2);
2609   ltrceC->fVecGZ=new TVectorD(557568/2);
2610   ltrceC->fVecLX=new TVectorD(557568/2);
2611   ltrceC->fVecLY=new TVectorD(557568/2);
2612   ltrceC->fVecLZ=new TVectorD(557568/2);
2613   arr->AddAt(ltrceC,337); //CE C-Side
2614
2615   //Calculate ideal positoins
2616   Float_t gxyz[3];
2617   Float_t lxyz[3];
2618   Int_t channelSideA=0;
2619   Int_t channelSideC=0;
2620   Int_t channelSide=0;
2621   AliTPCLaserTrack *ltrce=0x0;
2622
2623   for (Int_t isector=0; isector<72; ++isector){
2624     Int_t side=((isector/18)%2);
2625     for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2626       for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2627         fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2628         fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2629         if (side==0) {
2630           ltrce=ltrceA;
2631           channelSide=channelSideA;
2632         } else {
2633           ltrce=ltrceC;
2634           channelSide=channelSideC;
2635         }
2636
2637         ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2638         ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2639         ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2640         ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2641         ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2642 //         ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2643         ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2644         ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2645 //         ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2646
2647         if (side==0){
2648           ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2649           ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2650           ++channelSideA;
2651         }
2652         else {
2653           ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2654           ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2655           ++channelSideC;
2656         }
2657       }
2658     }
2659   }
2660
2661
2662 }
2663
2664 //_____________________________________________________________________
2665 void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2666 {
2667   //
2668   // calculate the drift velocity from the reconstructed clusters associated
2669   // to the ideal laser tracks
2670   // use two different fit scenarios: Separate fit for A- and C-Side
2671   //                                  Common fit for A- and C-Side
2672   //
2673
2674   if (!fArrFitGraphs){
2675     fArrFitGraphs=new TObjArray;
2676   }
2677
2678 //   static TLinearFitter fdriftA(5,"hyp4");
2679 //   static TLinearFitter fdriftC(5,"hyp4");
2680 //   static TLinearFitter fdriftAC(6,"hyp5");
2681   Double_t timestamp=fTimeBursts[burst];
2682
2683   static TLinearFitter fdriftA(4,"hyp3");
2684   static TLinearFitter fdriftC(4,"hyp3");
2685   static TLinearFitter fdriftAC(5,"hyp4");
2686   TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2687
2688   Float_t chi2A = 10;
2689   Float_t chi2C = 10;
2690   Float_t chi2AC = 10;
2691   Int_t npointsA=0;
2692   Int_t npointsC=0;
2693   Int_t npointsAC=0;
2694
2695   Double_t minres[3]={20.,1,0.8};
2696   //----
2697   for(Int_t i=0;i<3;i++){
2698
2699     fdriftA.ClearPoints();
2700     fdriftC.ClearPoints();
2701     fdriftAC.ClearPoints();
2702
2703     chi2A = 10;
2704     chi2C = 10;
2705     chi2AC = 10;
2706     npointsA=0;
2707     npointsC=0;
2708     npointsAC=0;
2709
2710     for (Int_t itrack=0; itrack<338; ++itrack){
2711       AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2712       AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2713
2714       //-- Exclude the tracks which has the biggest inclanation angle
2715       if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2716       Int_t clustercounter=0;
2717       Int_t indexMax=159;
2718
2719       //-- exclude the low intensity tracks
2720
2721       for (Int_t index=0; index<indexMax; ++index){
2722
2723         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2724         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2725         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2726
2727         if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2728       }
2729       if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2730       clustercounter=0;
2731
2732
2733       //-- drift length
2734       Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2735
2736       if (itrack>335) indexMax=557568/2;
2737       for (Int_t index=0; index<indexMax; ++index){
2738         Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2739         Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2740         Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2741         Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2742
2743         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2744         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2745         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2746         Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2747
2748       //cut if no track info available
2749         if (iltr->GetBundle()==0) continue;
2750         if (iR<133||mR<133) continue;
2751         if(mltr->fVecP2->GetMatrixArray()[index]>minres[i]) continue;
2752
2753         Double_t ldrift  = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2754         Double_t mdrift  = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2755
2756         //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2757         Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2758
2759         if (iltr->GetSide()==0){
2760           fdriftA.AddPoint(xxx,mdrift,1);
2761         }else{
2762           fdriftC.AddPoint(xxx,mdrift,1);
2763         }
2764 //         Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2765         Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2766         fdriftAC.AddPoint(xxx2,mdrift,1);
2767
2768       }//end index loop
2769     }//end laser track loop
2770
2771   //perform fit
2772     fdriftA.Eval();
2773     fdriftC.Eval();
2774     fdriftAC.Eval();
2775
2776
2777
2778   //get fit values
2779     fdriftA.GetParameters(fitA);
2780     fdriftC.GetParameters(fitC);
2781     fdriftAC.GetParameters(fitAC);
2782
2783   //Parameters:  0 linear offset
2784   //             1 mean drift velocity correction factor
2785   //             2 relative global y gradient
2786   //             3 radial deformation
2787   //             4 IROC/OROC offset
2788
2789 //      FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2790
2791     for (Int_t itrack=0; itrack<338; ++itrack){
2792       AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2793       AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2794
2795       //-- Exclude the tracks which has the biggest inclanation angle
2796       if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2797       Int_t clustercounter=0;
2798       Int_t indexMax=159;
2799
2800       //-- exclude the low intensity tracks
2801
2802       for (Int_t index=0; index<indexMax; ++index){
2803         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2804         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2805         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2806         if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2807       }
2808       if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2809       clustercounter=0;
2810
2811       //-- drift length
2812       Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2813
2814       if (itrack>335) indexMax=557568/2;
2815       for (Int_t index=0; index<indexMax; ++index){
2816         Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2817         Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2818         Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2819         Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2820
2821         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2822         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2823         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2824         Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2825
2826       //cut if no track info available
2827         if (iR<60||mR<60) continue;
2828
2829         Double_t ldrift  = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2830         Double_t mdrift  = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2831
2832         TVectorD *v=&fitA;
2833         if (iltr->GetSide()==1) v=&fitC;
2834 //         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);
2835         Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2836
2837         mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2838
2839       }
2840     }
2841
2842     fitA.ResizeTo(7);
2843     fitC.ResizeTo(7);
2844     fitAC.ResizeTo(8);
2845
2846 //set statistics values
2847
2848     npointsA= fdriftA.GetNpoints();
2849     if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2850     fitA[5]=npointsA;
2851     fitA[6]=chi2A;
2852
2853     npointsC= fdriftC.GetNpoints();
2854     if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2855     fitC[5]=npointsC;
2856     fitC[6]=chi2C;
2857
2858     npointsAC= fdriftAC.GetNpoints();
2859     if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2860     fitAC[5]=npointsAC;
2861     fitAC[6]=chi2AC;
2862
2863     if (fStreamLevel>2){
2864     //laser track information
2865       (*GetDebugStreamer()) << "DriftV" <<
2866         "iter="   << i <<
2867         "run="    << fRunNumber <<
2868         "time="   << timestamp <<
2869         "fitA.="  << &fitA <<
2870         "fitC.="  << &fitC <<
2871         "fitAC.=" << &fitAC <<
2872         "\n";
2873
2874
2875     }
2876
2877   }
2878 //-----
2879
2880
2881   //Parameters:  0 linear offset (global)
2882   //             1 mean drift velocity correction factor
2883   //             2 relative global y gradient
2884   //             3 radial deformation
2885   //             4 IROC/OROC offset
2886   //             5 linear offset relative A-C
2887
2888   //get graphs
2889   TGraphErrors *grA[7];
2890   TGraphErrors *grC[7];
2891   TGraphErrors *grAC[8];
2892   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_");
2893   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_");
2894
2895   TObjArray *arrNames=names.Tokenize(";");
2896   TObjArray *arrNamesAC=namesAC.Tokenize(";");
2897
2898   //A-Side graphs
2899   for (Int_t i=0; i<7; ++i){
2900     TString grName=arrNames->UncheckedAt(i)->GetName();
2901     grName+="A";
2902     grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2903     if (!grA[i]){
2904       grA[i]=new TGraphErrors;
2905       grA[i]->SetName(grName.Data());
2906       grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2907       fArrFitGraphs->Add(grA[i]);
2908     }
2909 //     Int_t ipoint=grA[i]->GetN();
2910     Int_t ipoint=burst;
2911     grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2912     grA[i]->SetPointError(ipoint,60,.0001);
2913     if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2914   }
2915
2916   //C-Side graphs
2917   for (Int_t i=0; i<7; ++i){
2918     TString grName=arrNames->UncheckedAt(i)->GetName();
2919     grName+="C";
2920     grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2921     if (!grC[i]){
2922       grC[i]=new TGraphErrors;
2923       grC[i]->SetName(grName.Data());
2924       grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2925       fArrFitGraphs->Add(grC[i]);
2926     }
2927 //     Int_t ipoint=grC[i]->GetN();
2928     Int_t ipoint=burst;
2929     grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2930     grC[i]->SetPointError(ipoint,60,.0001);
2931     if (i<4) grC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2932   }
2933
2934   //AC-Side graphs
2935   for (Int_t i=0; i<8; ++i){
2936     TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2937     grName+="AC";
2938     grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2939     if (!grAC[i]){
2940       grAC[i]=new TGraphErrors;
2941       grAC[i]->SetName(grName.Data());
2942       grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2943       fArrFitGraphs->Add(grAC[i]);
2944     }
2945 //     Int_t ipoint=grAC[i]->GetN();
2946     Int_t ipoint=burst;
2947     grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2948     grAC[i]->SetPointError(ipoint,60,.0001);
2949     if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2950   }
2951
2952   if (fDebugLevel>10){
2953     printf("A side fit parameters:\n");
2954     fitA.Print();
2955     printf("\nC side fit parameters:\n");
2956     fitC.Print();
2957     printf("\nAC side fit parameters:\n");
2958     fitAC.Print();
2959   }
2960   delete arrNames;
2961   delete arrNamesAC;
2962 }
2963
2964 //_____________________________________________________________________
2965 Double_t AliTPCCalibCE::SetBurstHnDrift()
2966 {
2967   //
2968   // Create a new THnSparse for the current burst
2969   // return the time of the current burst
2970   //
2971   Int_t i=0;
2972   for(i=0; i<fTimeBursts.GetNrows(); ++i){
2973     if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2974     if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2975       fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2976       return fTimeBursts(i);
2977     }
2978   }
2979
2980   CreateDVhist();
2981   fArrHnDrift.AddAt(fHnDrift,i);
2982   fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
2983 //   fPeaks[4]=0;
2984   return fTimeStamp;
2985 }
2986
2987 //_____________________________________________________________________
2988 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
2989 {
2990   //
2991   //  Write class to file
2992   //  option can be specified in the dir option:
2993   //  options:
2994   //    name=<objname>: the name of the calibration object in file will be <objname>
2995   //    type=<type>:    the saving type:
2996   //                    0 - write the complte object
2997   //                    1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
2998   //                    2 - like 2, but in addition delete objects that will most probably not be used for calibration
2999   //                    3 - store only calibration output, don't store the reference histograms
3000   //                        and THnSparse (requires Analyse called before)
3001   //
3002   //  NOTE: to read the object back, the ReadFromFile function should be used
3003   //
3004
3005   TString sDir(dir);
3006   TString objName=GetName();
3007   Int_t type=0;
3008
3009   //get options
3010   TObjArray *arr=sDir.Tokenize(",");
3011   TIter next(arr);
3012   TObjString *s;
3013   while ( (s=(TObjString*)next()) ){
3014     TString optString=s->GetString();
3015     optString.Remove(TString::kBoth,' ');
3016     if (optString.BeginsWith("name=")){
3017       objName=optString.ReplaceAll("name=","");
3018     }
3019     if (optString.BeginsWith("type=")){
3020       optString.ReplaceAll("type=","");
3021       type=optString.Atoi();
3022     }
3023   }
3024
3025   if (type==1||type==2) {
3026     //delete most probably not needed stuff
3027     //This requires Analyse to be called after reading the object from file
3028     fCalRocArrayT0.Delete();
3029     fCalRocArrayT0Err.Delete();
3030     fCalRocArrayQ.Delete();
3031     fCalRocArrayRMS.Delete();
3032     fCalRocArrayOutliers.Delete();
3033   }
3034   if (type==2||type==3){
3035     fParamArrayEventPol1.Delete();
3036     fParamArrayEventPol2.Delete();
3037   }
3038
3039   TObjArray histoQArray(72);
3040   TObjArray histoT0Array(72);
3041   TObjArray histoRMSArray(72);
3042   TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3043
3044   //save all large 2D histograms in separte pointers
3045   //to have a smaller memory print when saving the object
3046   if (type==1||type==2||type==3){
3047     for (Int_t i=0; i<72; ++i){
3048       histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3049       histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3050       histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3051     }
3052     fHistoQArray.SetOwner(kFALSE);
3053     fHistoT0Array.SetOwner(kFALSE);
3054     fHistoRMSArray.SetOwner(kFALSE);
3055     fHistoQArray.Clear();
3056     fHistoT0Array.Clear();
3057     fHistoRMSArray.Clear();
3058
3059     for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3060       arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3061     }
3062     fArrHnDrift.SetOwner(kFALSE);
3063     fArrHnDrift.Clear();
3064   }
3065
3066
3067   TDirectory *backup = gDirectory;
3068
3069   TFile f(filename,"recreate");
3070   this->Write(objName.Data());
3071   if (type==1||type==2) {
3072     histoQArray.Write("histoQArray",TObject::kSingleKey);
3073     histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3074     histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3075     arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3076   }
3077
3078   f.Save();
3079   f.Close();
3080
3081   //move histograms back to the object
3082   if (type==1||type==2){
3083     for (Int_t i=0; i<72; ++i){
3084       fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3085       fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3086       fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3087     }
3088     fHistoQArray.SetOwner(kTRUE);
3089     fHistoT0Array.SetOwner(kTRUE);
3090     fHistoRMSArray.SetOwner(kTRUE);
3091
3092     for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3093       fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3094     }
3095     fArrHnDrift.SetOwner(kTRUE);
3096   }
3097
3098   if ( backup ) backup->cd();
3099 }
3100 //_____________________________________________________________________
3101 AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3102 {
3103   //
3104   // Read object from file
3105   // Handle properly if the histogram arrays were stored separately
3106   // call Analyse to make sure to have the calibration relevant information in the object
3107   //
3108
3109   TFile f(filename);
3110   if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3111   TList *l=f.GetListOfKeys();
3112   TIter next(l);
3113   TKey *key=0x0;
3114   TObject *o=0x0;
3115
3116   AliTPCCalibCE *ce=0x0;
3117   TObjArray *histoQArray=0x0;
3118   TObjArray *histoT0Array=0x0;
3119   TObjArray *histoRMSArray=0x0;
3120   TObjArray *arrHnDrift=0x0;
3121
3122   while ( (key=(TKey*)next()) ){
3123     o=key->ReadObj();
3124     if ( o->IsA()==AliTPCCalibCE::Class() ){
3125       ce=(AliTPCCalibCE*)o;
3126     } else if ( o->IsA()==TObjArray::Class() ){
3127       TString name=key->GetName();
3128       if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3129       if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3130       if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3131       if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3132     }
3133   }
3134
3135   if (ce){
3136   //move histograms back to the object
3137     TH1* hist=0x0;
3138     if (histoQArray){
3139       for (Int_t i=0; i<72; ++i){
3140         hist=(TH1*)histoQArray->UncheckedAt(i);
3141         if (hist) hist->SetDirectory(0x0);
3142         ce->fHistoQArray.AddAt(hist,i);
3143       }
3144       ce->fHistoQArray.SetOwner(kTRUE);
3145     }
3146
3147     if (histoT0Array) {
3148       for (Int_t i=0; i<72; ++i){
3149         hist=(TH1*)histoT0Array->UncheckedAt(i);
3150         if (hist) hist->SetDirectory(0x0);
3151         ce->fHistoT0Array.AddAt(hist,i);
3152       }
3153       ce->fHistoT0Array.SetOwner(kTRUE);
3154     }
3155
3156     if (histoRMSArray){
3157       for (Int_t i=0; i<72; ++i){
3158         hist=(TH1*)histoRMSArray->UncheckedAt(i);
3159         if (hist) hist->SetDirectory(0x0);
3160         ce->fHistoRMSArray.AddAt(hist,i);
3161       }
3162       ce->fHistoRMSArray.SetOwner(kTRUE);
3163     }
3164
3165     if (arrHnDrift){
3166       for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3167         THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3168         ce->fArrHnDrift.AddAt(hSparse,i);
3169       }
3170     }
3171
3172     ce->Analyse();
3173   }
3174   f.Close();
3175
3176   return ce;
3177 }