1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////////////////
20 // Implementation of the TPC Central Electrode calibration //
22 // Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
24 ////////////////////////////////////////////////////////////////////////////////////////
27 // *************************************************************************************
28 // * Class Description *
29 // *************************************************************************************
32 <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
33 using laser runs.</h4>
35 The information retrieved is
36 <ul style="list-style-type: square;">
37 <li>Time arrival from the CE</li>
43 <ol style="list-style-type: upper-roman;">
44 <li><a href="#working">Working principle</a></li>
45 <li><a href="#user">User interface for filling data</a></li>
46 <li><a href="#info">Stored information</a></li>
49 <h3><a name="working">I. Working principle</a></h3>
51 <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
52 (see below). These in the end call the Update(...) function.</h4>
54 <ul style="list-style-type: square;">
55 <li>the Update(...) function:<br />
56 In this function the array fPadSignal is filled with the adc signals between the specified range
57 fFirstTimeBin and fLastTimeBin for the current pad.
58 before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
61 <ul style="list-style-type: square;">
62 <li>the ProcessPad() function:</li>
63 <ol style="list-style-type: decimal;">
64 <li>Find Pedestal and Noise information</li>
65 <ul style="list-style-type: square;">
66 <li>use database information which has to be set by calling<br />
67 SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
68 <li>if no information from the pedestal data base
69 is available the informaion is calculated on the fly
70 ( see FindPedestal() function )</li>
72 <li>Find local maxima of the pad signal</li>
73 <ul style="list-style-type: square;">
74 <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
75 have been observed ( see FindLocalMaxima(...) )</li>
77 <li>Find the CE signal information</li>
78 <ul style="list-style-type: square;">
79 <li>to find the position of the CE signal the Tmean information from the previos event is used
80 as the CE signal the local maximum closest to this Tmean is identified</li>
81 <li>calculate mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
82 the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
84 <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
85 <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
90 <h4>At the end of each event the EndEvent() function is called</h4>
92 <ul style="list-style-type: square;">
93 <li>the EndEvent() function:</li>
94 <ul style="list-style-type: square;">
95 <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
96 This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
97 <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
98 <li>calculate Mean Q</li>
99 <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
103 <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
104 <ul style="list-style-type: square;">
105 <li>the Analyse() function:</li>
106 <ul style="list-style-type: square;">
107 <li>calculate the mean values of T0, RMS, Q for each pad, using
108 the AliMathBase::GetCOG(...) function</li>
109 <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
110 (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
114 <h3><a name="user">II. User interface for filling data</a></h3>
116 <h4>To Fill information one of the following functions can be used:</h4>
118 <ul style="list-style-type: none;">
119 <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
120 <ul style="list-style-type: square;">
121 <li>process Date event</li>
122 <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
126 <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
127 <ul style="list-style-type: square;">
128 <li>process AliRawReader event</li>
129 <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
133 <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
134 <ul style="list-style-type: square;">
135 <li>process event from AliTPCRawStream</li>
136 <li>call Update function for signal filling</li>
140 <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
141 iPad, const Int_t iTimeBin, const Float_t signal);</li>
142 <ul style="list-style-type: square;">
143 <li>directly fill signal information (sector, row, pad, time bin, pad)
144 to the reference histograms</li>
148 <h4>It is also possible to merge two independently taken calibrations using the function</h4>
150 <ul style="list-style-type: none;">
151 <li> void Merge(AliTPCCalibSignal *sig)</li>
152 <ul style="list-style-type: square;">
153 <li>copy histograms in 'sig' if they do not exist in this instance</li>
154 <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
155 <li>After merging call Analyse again!</li>
160 <h4>example: filling data using root raw data:</h4>
162 void fillCE(Char_t *filename)
164 rawReader = new AliRawReaderRoot(fileName);
165 if ( !rawReader ) return;
166 AliTPCCalibCE *calib = new AliTPCCalibCE;
167 while (rawReader->NextEvent()){
168 calib->ProcessEvent(rawReader);
171 calib->DumpToFile("CEData.root");
177 <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
179 <h4><a name="info:stored">III.1 Stored information</a></h4>
180 <ul style="list-style-type: none;">
182 <ul style="list-style-type: none;">
183 <li>For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum])
184 is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
185 stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
189 <li>Calibration Data:</li>
190 <ul style="list-style-type: none;">
191 <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
192 the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
193 fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
194 is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
198 <li>For each event the following information is stored:</li>
200 <ul style="list-style-type: square;">
201 <li>event time ( TVectorD fVEventTime )</li>
202 <li>event id ( TVectorD fVEventNumber )</li>
204 <li>mean arrival time for each ROC ( TObjArray fTMeanArrayEvent )</li>
205 <li>mean Q for each ROC ( TObjArray fQMeanArrayEvent )</li>
206 <li>parameters of a plane fit for each ROC ( TObjArray fParamArrayEventPol1 )</li>
207 <li>parameters of a 2D parabola fit for each ROC ( TObjArray fParamArrayEventPol2 )</li>
211 <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
212 <ul style="list-style-type: none;">
213 <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
214 <ul style="list-style-type: square;">
215 <li>TH2F *GetHistoT0(Int_t sector);</li>
216 <li>TH2F *GetHistoRMS(Int_t sector);</li>
217 <li>TH2F *GetHistoQ(Int_t sector);</li>
221 <li>Accessing the calibration storage objects:</li>
222 <ul style="list-style-type: square;">
223 <li>AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values</li>
224 <li>AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values</li>
225 <li>AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values</li>
229 <li>Accessin the event by event information:</li>
230 <ul style="list-style-type: square;">
231 <li>The event by event information can be displayed using the</li>
232 <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
233 <li>which creates a graph from the specified variables</li>
237 <h4>example for visualisation:</h4>
239 //if the file "CEData.root" was created using the above example one could do the following:
240 TFile fileCE("CEData.root")
241 AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
242 ce->GetCalRocT0(0)->Draw("colz");
243 ce->GetCalRocRMS(0)->Draw("colz");
245 //or use the AliTPCCalPad functionality:
246 AliTPCCalPad padT0(ped->GetCalPadT0());
247 AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
248 padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
249 padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
251 //display event by event information:
252 //Draw mean arrival time as a function of the event time for oroc sector A00
253 ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
254 //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
255 ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");
258 //////////////////////////////////////////////////////////////////////////////////////
262 #include <TObjArray.h>
266 #include <TVectorF.h>
267 #include <TVectorD.h>
268 #include <TMatrixD.h>
273 #include <TDirectory.h>
276 #include <TCollection.h>
280 #include "AliRawReader.h"
281 #include "AliRawReaderRoot.h"
282 #include "AliRawReaderDate.h"
283 #include "AliRawEventHeaderBase.h"
284 #include "AliTPCRawStream.h"
285 #include "AliTPCRawStreamFast.h"
286 #include "AliTPCcalibDB.h"
287 #include "AliTPCCalROC.h"
288 #include "AliTPCCalPad.h"
289 #include "AliTPCROC.h"
290 #include "AliTPCParam.h"
291 #include "AliTPCCalibCE.h"
292 #include "AliMathBase.h"
293 #include "TTreeStream.h"
297 ClassImp(AliTPCCalibCE)
300 AliTPCCalibCE::AliTPCCalibCE() :
301 AliTPCCalibRawBase(),
315 fNoiseThresholdMax(5.),
316 fNoiseThresholdSum(8.),
317 fIsZeroSuppressed(kFALSE),
320 fParam(new AliTPCParam),
326 fCalRocArrayT0Err(72),
329 fCalRocArrayOutliers(72),
337 fParamArrayEventPol1(72),
338 fParamArrayEventPol2(72),
339 fTMeanArrayEvent(72),
340 fQMeanArrayEvent(72),
347 fPadTimesArrayEvent(72),
349 fPadRMSArrayEvent(72),
350 fPadPedestalArrayEvent(72),
360 fVTime0OffsetCounter(72),
366 // AliTPCSignal default constructor
368 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
372 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
374 //_____________________________________________________________________
375 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
376 AliTPCCalibRawBase(sig),
377 fNbinsT0(sig.fNbinsT0),
378 fXminT0(sig.fXminT0),
379 fXmaxT0(sig.fXmaxT0),
380 fNbinsQ(sig.fNbinsQ),
383 fNbinsRMS(sig.fNbinsRMS),
384 fXminRMS(sig.fXminRMS),
385 fXmaxRMS(sig.fXmaxRMS),
386 fPeakDetMinus(sig.fPeakDetMinus),
387 fPeakDetPlus(sig.fPeakDetPlus),
388 fPeakIntMinus(sig.fPeakIntMinus),
389 fPeakIntPlus(sig.fPeakIntPlus),
390 fNoiseThresholdMax(sig.fNoiseThresholdMax),
391 fNoiseThresholdSum(sig.fNoiseThresholdSum),
392 fIsZeroSuppressed(sig.fIsZeroSuppressed),
395 fParam(new AliTPCParam),
401 fCalRocArrayT0Err(72),
404 fCalRocArrayOutliers(72),
408 fMeanT0rms(sig.fMeanT0rms),
409 fMeanQrms(sig.fMeanQrms),
410 fMeanRMSrms(sig.fMeanRMSrms),
412 fParamArrayEventPol1(72),
413 fParamArrayEventPol2(72),
414 fTMeanArrayEvent(72),
415 fQMeanArrayEvent(72),
416 fVEventTime(sig.fVEventTime),
417 fVEventNumber(sig.fVEventNumber),
418 fVTime0SideA(sig.fVTime0SideA),
419 fVTime0SideC(sig.fVTime0SideC),
422 fPadTimesArrayEvent(72),
424 fPadRMSArrayEvent(72),
425 fPadPedestalArrayEvent(72),
435 fVTime0OffsetCounter(72),
441 // AliTPCSignal copy constructor
443 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
445 for (Int_t iSec = 0; iSec < 72; ++iSec){
446 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
447 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
448 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
449 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
451 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
452 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
453 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
455 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
456 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
457 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
458 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
461 TH2S *hNew = new TH2S(*hQ);
462 hNew->SetDirectory(0);
463 fHistoQArray.AddAt(hNew,iSec);
466 TH2S *hNew = new TH2S(*hT0);
467 hNew->SetDirectory(0);
468 fHistoT0Array.AddAt(hNew,iSec);
471 TH2S *hNew = new TH2S(*hRMS);
472 hNew->SetDirectory(0);
473 fHistoRMSArray.AddAt(hNew,iSec);
477 //copy fit parameters event by event
479 for (Int_t iSec=0; iSec<72; ++iSec){
480 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
482 TObjArray *arrEvents = new TObjArray(arr->GetSize());
483 fParamArrayEventPol1.AddAt(arrEvents, iSec);
484 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
485 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
486 arrEvents->AddAt(new TVectorD(*vec),iEvent);
489 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
491 TObjArray *arrEvents = new TObjArray(arr->GetSize());
492 fParamArrayEventPol2.AddAt(arrEvents, iSec);
493 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
494 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
495 arrEvents->AddAt(new TVectorD(*vec),iEvent);
498 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
499 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
501 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
503 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
507 fVEventTime.ResizeTo(sig.fVEventTime);
508 fVEventNumber.ResizeTo(sig.fVEventNumber);
509 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
510 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
514 //_____________________________________________________________________
515 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
516 AliTPCCalibRawBase(),
530 fNoiseThresholdMax(5.),
531 fNoiseThresholdSum(8.),
532 fIsZeroSuppressed(kFALSE),
535 fParam(new AliTPCParam),
541 fCalRocArrayT0Err(72),
544 fCalRocArrayOutliers(72),
552 fParamArrayEventPol1(72),
553 fParamArrayEventPol2(72),
554 fTMeanArrayEvent(72),
555 fQMeanArrayEvent(72),
562 fPadTimesArrayEvent(72),
564 fPadRMSArrayEvent(72),
565 fPadPedestalArrayEvent(72),
575 fVTime0OffsetCounter(72),
581 // constructor which uses a tmap as input to set some specific parameters
583 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
586 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
587 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
588 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
589 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
590 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
591 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
592 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
593 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
594 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
595 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
596 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
597 if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
598 if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
599 if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
600 if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
601 if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
602 if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
603 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
604 if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
605 if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
607 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
612 //_____________________________________________________________________
613 AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
616 // assignment operator
618 if (&source == this) return *this;
619 new (this) AliTPCCalibCE(source);
623 //_____________________________________________________________________
624 AliTPCCalibCE::~AliTPCCalibCE()
630 fCalRocArrayT0.Delete();
631 fCalRocArrayT0Err.Delete();
632 fCalRocArrayQ.Delete();
633 fCalRocArrayRMS.Delete();
634 fCalRocArrayOutliers.Delete();
636 fHistoQArray.Delete();
637 fHistoT0Array.Delete();
638 fHistoRMSArray.Delete();
640 fHistoTmean.Delete();
642 fParamArrayEventPol1.Delete();
643 fParamArrayEventPol2.Delete();
644 fTMeanArrayEvent.Delete();
645 fQMeanArrayEvent.Delete();
647 fPadTimesArrayEvent.Delete();
648 fPadQArrayEvent.Delete();
649 fPadRMSArrayEvent.Delete();
650 fPadPedestalArrayEvent.Delete();
652 // if ( fHTime0 ) delete fHTime0;
655 //_____________________________________________________________________
656 Int_t AliTPCCalibCE::Update(const Int_t icsector,
659 const Int_t icTimeBin,
660 const Float_t csignal)
663 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
664 // no extra analysis necessary. Assumes knowledge of the signal shape!
665 // assumes that it is looped over consecutive time bins of one pad
670 if (icRow<0) return 0;
671 if (icPad<0) return 0;
672 if (icTimeBin<0) return 0;
673 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
675 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
677 //init first pad and sector in this event
678 if ( fCurrentChannel == -1 ) {
680 fCurrentChannel = iChannel;
681 fCurrentSector = icsector;
685 //process last pad if we change to a new one
686 if ( iChannel != fCurrentChannel ){
688 fLastSector=fCurrentSector;
689 fCurrentChannel = iChannel;
690 fCurrentSector = icsector;
694 //fill signals for current pad
695 fPadSignal[icTimeBin]=csignal;
696 if ( csignal > fMaxPadSignal ){
697 fMaxPadSignal = csignal;
698 fMaxTimeBin = icTimeBin;
702 //_____________________________________________________________________
703 void AliTPCCalibCE::FindPedestal(Float_t part)
706 // find pedestal and noise for the current pad. Use either database or
707 // truncated mean with part*100%
709 Bool_t noPedestal = kTRUE;
711 //use pedestal database if set
712 if (fPedestalTPC&&fPadNoiseTPC){
713 //only load new pedestals if the sector has changed
714 if ( fCurrentSector!=fLastSector ){
715 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
716 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
719 if ( fPedestalROC&&fPadNoiseROC ){
720 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
721 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
727 //if we are not running with pedestal database, or for the current sector there is no information
728 //available, calculate the pedestal and noise on the fly
732 if ( fIsZeroSuppressed ) return;
733 const Int_t kPedMax = 100; //maximum pedestal value
742 UShort_t histo[kPedMax];
743 memset(histo,0,kPedMax*sizeof(UShort_t));
745 //fill pedestal histogram
746 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
747 padSignal = fPadSignal[i];
748 if (padSignal<=0) continue;
749 if (padSignal>max && i>10) {
753 if (padSignal>kPedMax-1) continue;
754 histo[int(padSignal+0.5)]++;
758 for (Int_t i=1; i<kPedMax; ++i){
759 if (count1<count0*0.5) median=i;
764 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
766 for (Int_t idelta=1; idelta<10; ++idelta){
767 if (median-idelta<=0) continue;
768 if (median+idelta>kPedMax) continue;
769 if (count<part*count1){
770 count+=histo[median-idelta];
771 mean +=histo[median-idelta]*(median-idelta);
772 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
773 count+=histo[median+idelta];
774 mean +=histo[median+idelta]*(median+idelta);
775 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
780 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
786 //_____________________________________________________________________
787 void AliTPCCalibCE::UpdateCETimeRef()
789 // Find the time reference of the last valid CE signal in sector
790 // for irocs of the A-Side the reference of the corresponging OROC is returned
791 // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
792 if ( fLastSector == fCurrentSector ) return;
793 Int_t sector=fCurrentSector;
794 if ( sector < 18 ) sector+=36;
796 TVectorF *vtRef = GetTMeanEvents(sector);
797 if ( !vtRef ) return;
798 Int_t vtRefSize= vtRef->GetNrows();
799 if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
800 else vtRefSize=fNevents;
801 while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
802 fCurrentCETimeRef=(*vtRef)[vtRefSize];
803 AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
805 //_____________________________________________________________________
806 void AliTPCCalibCE::FindCESignal(TVectorD ¶m, Float_t &qSum, const TVectorF maxima)
809 // Find position, signal width and height of the CE signal (last signal)
810 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
811 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
814 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
816 Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise; // threshold for the signal sum
817 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
818 const Int_t kCemax = fPeakIntPlus;
820 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
822 // find maximum closest to the sector mean from the last event
823 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
824 // get sector mean of last event
825 Float_t tmean = fCurrentCETimeRef;
826 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
827 minDist = tmean-maxima[imax];
828 cemaxpos = (Int_t)maxima[imax];
831 // printf("L1 phase TB: %f\n",GetL1PhaseTB());
833 ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
834 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
835 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
836 Float_t signal = fPadSignal[i]-fPadPedestal;
838 ceTime+=signal*(i+0.5);
839 ceRMS +=signal*(i+0.5)*(i+0.5);
845 if (ceQmax&&ceQsum>ceSumThreshold) {
847 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
848 ceTime-=GetL1PhaseTB();
849 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
850 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
852 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
853 // the pick-up signal should scale with the pad area. In addition
854 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
855 // ratio 2/3. The pad area we express in cm2. We normalise the signal
856 // to the OROC signal (factor 2/3 for the IROCs).
857 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
858 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
861 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
862 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
874 //_____________________________________________________________________
875 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
878 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
879 // and 'tplus' timebins after 'pos'
881 if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
882 for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
883 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
884 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
885 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
887 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
891 //_____________________________________________________________________
892 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
895 // Find local maxima on the pad signal and Histogram them
897 Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal
900 for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
901 if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
902 if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
903 if (count<maxima.GetNrows()){
904 maxima.GetMatrixArray()[count++]=i;
905 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
906 i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin fPeakDetMinus+fPeakDetPlus-1
911 //_____________________________________________________________________
912 void AliTPCCalibCE::ProcessPad()
915 // Process data of current pad
919 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
920 // + central electrode and possibly post peaks from the CE signal
921 // however if we are on a high noise pad a lot more peaks due to the noise might occur
922 FindLocalMaxima(maxima);
923 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
925 UpdateCETimeRef(); // update the time refenrence for the current sector
926 if ( fCurrentCETimeRef<1e-30 ) return; //return if we don't have time 0 info, eg if only one side has laser
929 FindCESignal(param, qSum, maxima);
931 Double_t meanT = param[1];
932 Double_t sigmaT = param[2];
934 //Fill Event T0 counter
935 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
938 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
941 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
944 //Fill debugging info
945 if ( GetStreamLevel()>0 ){
946 (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
947 (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
948 (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
953 //_____________________________________________________________________
954 void AliTPCCalibCE::EndEvent()
956 // Process data of current pad
957 // The Functions 'SetTimeStamp' and 'SetRunNumber' should be called
958 // before the EndEvent function to set the event timestamp and number!!!
959 // This is automatically done if the ProcessEvent(AliRawReader *rawReader)
960 // function was called
962 //check if last pad has allready been processed, if not do so
963 if ( fMaxTimeBin>-1 ) ProcessPad();
965 AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
969 // TVectorF vMeanTime(72);
970 // TVectorF vMeanQ(72);
971 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
972 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
974 //find mean time0 offset for side A and C
975 //use only orocs due to the better statistics
976 Double_t time0Side[2]; //time0 for side A:0 and C:1
977 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:1
978 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
979 for ( Int_t iSec = 36; iSec<72; ++iSec ){
980 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
981 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
983 if ( time0SideCount[0] >0 )
984 time0Side[0]/=time0SideCount[0];
985 if ( time0SideCount[1] >0 )
986 time0Side[1]/=time0SideCount[1];
987 // end find time0 offset
988 AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
990 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
991 for ( Int_t iSec = 0; iSec<72; ++iSec ){
992 AliDebug(4,Form("Processing sector '%02d'\n",iSec));
993 //find median and then calculate the mean around it
994 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
995 if ( !hMeanT ) continue;
996 //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
997 if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
999 AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1003 Double_t entries = hMeanT->GetEffectiveEntries();
1005 Short_t *arr = hMeanT->GetArray()+1;
1007 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1009 if ( sum>=(entries/2.) ) break;
1012 Int_t firstBin = fFirstTimeBin+ibin-delta;
1013 Int_t lastBin = fFirstTimeBin+ibin+delta;
1014 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1015 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
1016 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1018 // check boundaries for ebye info of mean time
1019 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1020 Int_t vSize=vMeanTime->GetNrows();
1021 if ( vSize < fNevents+1 ){
1022 vMeanTime->ResizeTo(vSize+100);
1025 // store mean time for the readout sides
1026 vSize=fVTime0SideA.GetNrows();
1027 if ( vSize < fNevents+1 ){
1028 fVTime0SideA.ResizeTo(vSize+100);
1029 fVTime0SideC.ResizeTo(vSize+100);
1031 fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1032 fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
1034 vMeanTime->GetMatrixArray()[fNevents]=median;
1038 TVectorF *vTimes = GetPadTimesEvent(iSec);
1039 if ( !vTimes ) continue; //continue if no time information for this sector is available
1041 AliTPCCalROC calIrocOutliers(0);
1042 AliTPCCalROC calOrocOutliers(36);
1044 // calculate mean Q of the sector
1045 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1046 vSize=vMeanQ->GetNrows();
1047 if ( vSize < fNevents+1 ){
1048 vMeanQ->ResizeTo(vSize+100);
1051 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1052 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1054 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1055 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
1057 //set values for temporary roc calibration class
1059 calIroc->SetValue(iChannel, time);
1060 if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
1063 calOroc->SetValue(iChannel, time);
1064 if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
1067 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1068 // test that we really found the CE signal reliably
1069 if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1070 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
1074 //------------------------------- Debug start ------------------------------
1075 if ( GetStreamLevel()>0 ){
1076 TTreeSRedirector *streamer=GetDebugStreamer();
1082 Float_t q = (*GetPadQEvent(iSec))[iChannel];
1083 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1085 UInt_t channel=iChannel;
1088 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1089 pad = channel-fROC->GetRowIndexes(sector)[row];
1090 padc = pad-(fROC->GetNPads(sector,row)/2);
1092 // TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1093 // Form("hSignalD%d.%d.%d",sector,row,pad),
1094 // fLastTimeBin-fFirstTimeBin,
1095 // fFirstTimeBin,fLastTimeBin);
1096 // h1->SetDirectory(0);
1098 // for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1099 // h1->Fill(i,fPadSignal(i));
1102 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1103 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1104 Double_t t0Side = time0Side[(iSec/18)%2];
1105 (*streamer) << "DataPad" <<
1106 "Event=" << fNevents <<
1107 "RunNumber=" << fRunNumber <<
1108 "TimeStamp=" << fTimeStamp <<
1109 "Sector="<< sector <<
1113 "PadSec="<< channel <<
1114 "Time0Sec=" << t0Sec <<
1115 "Time0Side=" << t0Side <<
1119 "MeanQ=" << meanQ <<
1120 // "hist.=" << h1 <<
1126 //----------------------------- Debug end ------------------------------
1127 }// end channel loop
1130 //do fitting now only in debug mode
1131 if (GetDebugLevel()>0){
1132 TVectorD paramPol1(3);
1133 TVectorD paramPol2(6);
1134 TMatrixD matPol1(3,3);
1135 TMatrixD matPol2(6,6);
1139 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1141 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1142 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1144 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1145 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1148 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1149 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1152 //------------------------------- Debug start ------------------------------
1153 if ( GetStreamLevel()>0 ){
1154 TTreeSRedirector *streamer=GetDebugStreamer();
1156 (*streamer) << "DataRoc" <<
1157 // "Event=" << fEvent <<
1158 "RunNumber=" << fRunNumber <<
1159 "TimeStamp=" << fTimeStamp <<
1161 "hMeanT.=" << hMeanT <<
1162 "median=" << median <<
1163 "paramPol1.=" << ¶mPol1 <<
1164 "paramPol2.=" << ¶mPol2 <<
1165 "matPol1.=" << &matPol1 <<
1166 "matPol2.=" << &matPol2 <<
1167 "chi2Pol1=" << chi2Pol1 <<
1168 "chi2Pol2=" << chi2Pol2 <<
1173 //------------------------------- Debug end ------------------------------
1176 //return if no sector has a valid mean time
1177 if ( nSecMeanT == 0 ) return;
1180 // fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1181 // fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1182 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1183 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1184 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1186 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1187 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1190 fOldRunNumber = fRunNumber;
1194 AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1196 //_____________________________________________________________________
1197 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1198 Int_t nbinsY, Float_t ymin, Float_t ymax,
1199 const Char_t *type, Bool_t force)
1202 // return pointer to TH2S histogram of 'type'
1203 // if force is true create a new histogram if it doesn't exist allready
1205 if ( !force || arr->UncheckedAt(sector) )
1206 return (TH2S*)arr->UncheckedAt(sector);
1208 // if we are forced and histogram doesn't exist yet create it
1209 Char_t name[255], title[255];
1211 sprintf(name,"hCalib%s%.2d",type,sector);
1212 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1214 // new histogram with Q calib information. One value for each pad!
1215 TH2S* hist = new TH2S(name,title,
1217 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1218 hist->SetDirectory(0);
1219 arr->AddAt(hist,sector);
1222 //_____________________________________________________________________
1223 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1226 // return pointer to T0 histogram
1227 // if force is true create a new histogram if it doesn't exist allready
1229 TObjArray *arr = &fHistoT0Array;
1230 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1232 //_____________________________________________________________________
1233 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1236 // return pointer to Q histogram
1237 // if force is true create a new histogram if it doesn't exist allready
1239 TObjArray *arr = &fHistoQArray;
1240 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1242 //_____________________________________________________________________
1243 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1246 // return pointer to Q histogram
1247 // if force is true create a new histogram if it doesn't exist allready
1249 TObjArray *arr = &fHistoRMSArray;
1250 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1252 //_____________________________________________________________________
1253 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1254 const Char_t *type, Bool_t force)
1257 // return pointer to TH1S histogram
1258 // if force is true create a new histogram if it doesn't exist allready
1260 if ( !force || arr->UncheckedAt(sector) )
1261 return (TH1S*)arr->UncheckedAt(sector);
1263 // if we are forced and histogram doesn't yes exist create it
1264 Char_t name[255], title[255];
1266 sprintf(name,"hCalib%s%.2d",type,sector);
1267 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1269 // new histogram with calib information. One value for each pad!
1270 TH1S* hist = new TH1S(name,title,
1271 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1272 hist->SetDirectory(0);
1273 arr->AddAt(hist,sector);
1276 //_____________________________________________________________________
1277 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1280 // return pointer to Q histogram
1281 // if force is true create a new histogram if it doesn't exist allready
1283 TObjArray *arr = &fHistoTmean;
1284 return GetHisto(sector, arr, "LastTmean", force);
1286 //_____________________________________________________________________
1287 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1290 // return pointer to Pad Info from 'arr' for the current event and sector
1291 // if force is true create it if it doesn't exist allready
1293 if ( !force || arr->UncheckedAt(sector) )
1294 return (TVectorF*)arr->UncheckedAt(sector);
1296 TVectorF *vect = new TVectorF(size);
1297 arr->AddAt(vect,sector);
1300 //_____________________________________________________________________
1301 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1304 // return pointer to Pad Times Array for the current event and sector
1305 // if force is true create it if it doesn't exist allready
1307 TObjArray *arr = &fPadTimesArrayEvent;
1308 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1310 //_____________________________________________________________________
1311 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1314 // return pointer to Pad Q Array for the current event and sector
1315 // if force is true create it if it doesn't exist allready
1316 // for debugging purposes only
1319 TObjArray *arr = &fPadQArrayEvent;
1320 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1322 //_____________________________________________________________________
1323 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1326 // return pointer to Pad RMS Array for the current event and sector
1327 // if force is true create it if it doesn't exist allready
1328 // for debugging purposes only
1330 TObjArray *arr = &fPadRMSArrayEvent;
1331 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1333 //_____________________________________________________________________
1334 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1337 // return pointer to Pad RMS Array for the current event and sector
1338 // if force is true create it if it doesn't exist allready
1339 // for debugging purposes only
1341 TObjArray *arr = &fPadPedestalArrayEvent;
1342 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1344 //_____________________________________________________________________
1345 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1348 // return pointer to the EbyE info of the mean arrival time for 'sector'
1349 // if force is true create it if it doesn't exist allready
1351 TObjArray *arr = &fTMeanArrayEvent;
1352 return GetVectSector(sector,arr,100,force);
1354 //_____________________________________________________________________
1355 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1358 // return pointer to the EbyE info of the mean arrival time for 'sector'
1359 // if force is true create it if it doesn't exist allready
1361 TObjArray *arr = &fQMeanArrayEvent;
1362 return GetVectSector(sector,arr,100,force);
1364 //_____________________________________________________________________
1365 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1368 // return pointer to ROC Calibration
1369 // if force is true create a new histogram if it doesn't exist allready
1371 if ( !force || arr->UncheckedAt(sector) )
1372 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1374 // if we are forced and histogram doesn't yes exist create it
1376 // new AliTPCCalROC for T0 information. One value for each pad!
1377 AliTPCCalROC *croc = new AliTPCCalROC(sector);
1378 arr->AddAt(croc,sector);
1381 //_____________________________________________________________________
1382 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1385 // return pointer to Time 0 ROC Calibration
1386 // if force is true create a new histogram if it doesn't exist allready
1388 TObjArray *arr = &fCalRocArrayT0;
1389 return GetCalRoc(sector, arr, force);
1391 //_____________________________________________________________________
1392 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1395 // return pointer to the error of Time 0 ROC Calibration
1396 // if force is true create a new histogram if it doesn't exist allready
1398 TObjArray *arr = &fCalRocArrayT0Err;
1399 return GetCalRoc(sector, arr, force);
1401 //_____________________________________________________________________
1402 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1405 // return pointer to T0 ROC Calibration
1406 // if force is true create a new histogram if it doesn't exist allready
1408 TObjArray *arr = &fCalRocArrayQ;
1409 return GetCalRoc(sector, arr, force);
1411 //_____________________________________________________________________
1412 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1415 // return pointer to signal width ROC Calibration
1416 // if force is true create a new histogram if it doesn't exist allready
1418 TObjArray *arr = &fCalRocArrayRMS;
1419 return GetCalRoc(sector, arr, force);
1421 //_____________________________________________________________________
1422 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1425 // return pointer to Outliers
1426 // if force is true create a new histogram if it doesn't exist allready
1428 TObjArray *arr = &fCalRocArrayOutliers;
1429 return GetCalRoc(sector, arr, force);
1431 //_____________________________________________________________________
1432 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1435 // return pointer to TObjArray of fit parameters
1436 // if force is true create a new histogram if it doesn't exist allready
1438 if ( !force || arr->UncheckedAt(sector) )
1439 return (TObjArray*)arr->UncheckedAt(sector);
1441 // if we are forced and array doesn't yes exist create it
1443 // new TObjArray for parameters
1444 TObjArray *newArr = new TObjArray;
1445 arr->AddAt(newArr,sector);
1448 //_____________________________________________________________________
1449 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1452 // return pointer to TObjArray of fit parameters from plane fit
1453 // if force is true create a new histogram if it doesn't exist allready
1455 TObjArray *arr = &fParamArrayEventPol1;
1456 return GetParamArray(sector, arr, force);
1458 //_____________________________________________________________________
1459 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1462 // return pointer to TObjArray of fit parameters from parabola fit
1463 // if force is true create a new histogram if it doesn't exist allready
1465 TObjArray *arr = &fParamArrayEventPol2;
1466 return GetParamArray(sector, arr, force);
1468 //_____________________________________________________________________
1469 void AliTPCCalibCE::ResetEvent()
1472 // Reset global counters -- Should be called before each event is processed
1481 fPadTimesArrayEvent.Delete();
1482 fPadQArrayEvent.Delete();
1483 fPadRMSArrayEvent.Delete();
1484 fPadPedestalArrayEvent.Delete();
1486 for ( Int_t i=0; i<72; ++i ){
1487 fVTime0Offset.GetMatrixArray()[i]=0;
1488 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1489 fVMeanQ.GetMatrixArray()[i]=0;
1490 fVMeanQCounter.GetMatrixArray()[i]=0;
1493 //_____________________________________________________________________
1494 void AliTPCCalibCE::ResetPad()
1497 // Reset pad infos -- Should be called after a pad has been processed
1499 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1506 //_____________________________________________________________________
1507 void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
1510 // Merge ce to the current AliTPCCalibCE
1514 for (Int_t iSec=0; iSec<72; ++iSec){
1515 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1516 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1517 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1521 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1522 TH2S *hRefQ = GetHistoQ(iSec);
1523 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1525 TH2S *hist = new TH2S(*hRefQmerge);
1526 hist->SetDirectory(0);
1527 fHistoQArray.AddAt(hist, iSec);
1529 hRefQmerge->SetDirectory(dir);
1532 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1533 TH2S *hRefT0 = GetHistoT0(iSec);
1534 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1536 TH2S *hist = new TH2S(*hRefT0merge);
1537 hist->SetDirectory(0);
1538 fHistoT0Array.AddAt(hist, iSec);
1540 hRefT0merge->SetDirectory(dir);
1542 if ( hRefRMSmerge ){
1543 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1544 TH2S *hRefRMS = GetHistoRMS(iSec);
1545 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1547 TH2S *hist = new TH2S(*hRefRMSmerge);
1548 hist->SetDirectory(0);
1549 fHistoRMSArray.AddAt(hist, iSec);
1551 hRefRMSmerge->SetDirectory(dir);
1556 // merge time information
1559 Int_t nCEevents = ce->GetNeventsProcessed();
1560 for (Int_t iSec=0; iSec<72; ++iSec){
1561 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1562 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1563 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1564 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1566 TObjArray *arrPol1 = 0x0;
1567 TObjArray *arrPol2 = 0x0;
1568 TVectorF *vMeanTime = 0x0;
1569 TVectorF *vMeanQ = 0x0;
1572 if ( arrPol1CE && arrPol2CE ){
1573 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1574 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1575 arrPol1->Expand(fNevents+nCEevents);
1576 arrPol2->Expand(fNevents+nCEevents);
1578 if ( vMeanTimeCE && vMeanQCE ){
1579 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1580 vMeanQ = GetQMeanEvents(iSec,kTRUE);
1581 vMeanTime->ResizeTo(fNevents+nCEevents);
1582 vMeanQ->ResizeTo(fNevents+nCEevents);
1586 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1587 if ( arrPol1CE && arrPol2CE ){
1588 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1589 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1590 if ( paramPol1 && paramPol2 ){
1591 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1592 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1595 if ( vMeanTimeCE && vMeanQCE ){
1596 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1597 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1604 const TVectorD& eventTimes = ce->fVEventTime;
1605 const TVectorD& eventIds = ce->fVEventNumber;
1606 const TVectorF& time0SideA = ce->fVTime0SideA;
1607 const TVectorF& time0SideC = ce->fVTime0SideC;
1608 fVEventTime.ResizeTo(fNevents+nCEevents);
1609 fVEventNumber.ResizeTo(fNevents+nCEevents);
1610 fVTime0SideA.ResizeTo(fNevents+nCEevents);
1611 fVTime0SideC.ResizeTo(fNevents+nCEevents);
1613 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1614 Double_t evTime = eventTimes.GetMatrixArray()[iEvent];
1615 Double_t evId = eventIds.GetMatrixArray()[iEvent];
1616 Float_t t0SideA = time0SideA.GetMatrixArray()[iEvent];
1617 Float_t t0SideC = time0SideC.GetMatrixArray()[iEvent];
1619 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1620 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1621 fVTime0SideA.GetMatrixArray()[fNevents+iEvent] = t0SideA;
1622 fVTime0SideC.GetMatrixArray()[fNevents+iEvent] = t0SideC;
1624 fNevents+=nCEevents; //increase event counter
1627 //_____________________________________________________________________
1628 Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1631 // Merge all objects of this type in list
1637 AliTPCCalibCE *ce=0;
1640 while ( (o=next()) ){
1641 ce=dynamic_cast<AliTPCCalibCE*>(o);
1651 //_____________________________________________________________________
1652 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1655 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1656 // or side (-1: A-Side, -2: C-Side)
1657 // xVariable: 0-event time, 1-event id, 2-internal event counter
1658 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1659 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1660 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1661 // not used for mean time and mean Q )
1662 // for an example see class description at the beginning
1665 Double_t *x = new Double_t[fNevents];
1666 Double_t *y = new Double_t[fNevents];
1668 TVectorD *xVar = 0x0;
1669 TObjArray *aType = 0x0;
1673 if ( (sector<-2) || (sector>71) ) return 0x0; //sector outside valid range
1674 if ( (xVariable<0) || (xVariable>2) ) return 0x0; //invalid x-variable
1675 if ( (fitType<0) || (fitType>3) ) return 0x0; //invalid fit type
1676 if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1677 if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1678 if ( sector<0 && fitType!=2) return 0x0; //for side wise information only mean time is available
1682 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1683 aType = &fParamArrayEventPol1;
1684 if ( aType->At(sector)==0x0 ) return 0x0;
1686 else if ( fitType==1 ){
1687 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1688 aType = &fParamArrayEventPol2;
1689 if ( aType->At(sector)==0x0 ) return 0x0;
1693 if ( xVariable == 0 ) xVar = &fVEventTime;
1694 if ( xVariable == 1 ) xVar = &fVEventNumber;
1695 if ( xVariable == 2 ) {
1696 xVar = new TVectorD(fNevents);
1697 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1700 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1702 TObjArray *events = (TObjArray*)(aType->At(sector));
1703 if ( events->GetSize()<=ievent ) break;
1704 TVectorD *v = (TVectorD*)(events->At(ievent));
1705 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1706 } else if (fitType == 2) {
1707 Double_t xValue=(*xVar)[ievent];
1709 if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1710 else if (sector==-1) yValue=fVTime0SideA(ievent);
1711 else if (sector==-2) yValue=fVTime0SideC(ievent);
1712 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1713 }else if (fitType == 3) {
1714 Double_t xValue=(*xVar)[ievent];
1715 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1716 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1720 TGraph *gr = new TGraph(npoints);
1721 //sort xVariable increasing
1722 Int_t *sortIndex = new Int_t[npoints];
1723 TMath::Sort(npoints,x,sortIndex);
1724 for (Int_t i=0;i<npoints;++i){
1725 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1729 if ( xVariable == 2 ) delete xVar;
1732 delete [] sortIndex;
1735 //_____________________________________________________________________
1736 void AliTPCCalibCE::Analyse()
1739 // Calculate calibration constants
1743 TVectorD paramT0(3);
1744 TVectorD paramRMS(3);
1745 TMatrixD dummy(3,3);
1747 Float_t channelCounter=0;
1752 for (Int_t iSec=0; iSec<72; ++iSec){
1753 TH2S *hT0 = GetHistoT0(iSec);
1754 if (!hT0 ) continue;
1756 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1757 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1758 AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1759 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1760 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1762 TH2S *hQ = GetHistoQ(iSec);
1763 TH2S *hRMS = GetHistoRMS(iSec);
1765 Short_t *arrayhQ = hQ->GetArray();
1766 Short_t *arrayhT0 = hT0->GetArray();
1767 Short_t *arrayhRMS = hRMS->GetArray();
1769 UInt_t nChannels = fROC->GetNChannels(iSec);
1777 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1780 Float_t cogTime0 = -1000;
1781 Float_t cogQ = -1000;
1782 Float_t cogRMS = -1000;
1788 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1789 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1790 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1792 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
1794 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
1796 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
1801 //outlier specifications
1802 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1809 rocQ->SetValue(iChannel, cogQ*cogQ);
1810 rocT0->SetValue(iChannel, cogTime0);
1811 rocT0Err->SetValue(iChannel, rmsT0);
1812 rocRMS->SetValue(iChannel, cogRMS);
1813 rocOut->SetValue(iChannel, cogOut);
1817 if ( GetStreamLevel() > 0 ){
1818 TTreeSRedirector *streamer=GetDebugStreamer();
1821 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1822 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1823 padc = pad-(fROC->GetNPads(iSec,row)/2);
1825 (*streamer) << "DataEnd" <<
1826 "Sector=" << iSec <<
1830 "PadSec=" << iChannel <<
1832 "T0=" << cogTime0 <<
1842 if ( channelCounter>0 ){
1843 fMeanT0rms/=channelCounter;
1844 fMeanQrms/=channelCounter;
1845 fMeanRMSrms/=channelCounter;
1847 // if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
1848 // delete fDebugStreamer;
1849 // fDebugStreamer = 0x0;
1850 fVEventTime.ResizeTo(fNevents);
1851 fVEventNumber.ResizeTo(fNevents);
1852 fVTime0SideA.ResizeTo(fNevents);
1853 fVTime0SideC.ResizeTo(fNevents);