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