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