Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /// \class AliTPCCalibPulser
17 : /// \brief Implementation of the TPC pulser calibration
18 : ///
19 : /// \author Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
20 : ///
21 : /// Class Description
22 : /// -----------------
23 : ///
24 : /// The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
25 : /// runs performed with the calibration pulser.
26 : ///
27 : /// The information retrieved is
28 : /// - Time0 differences
29 : /// - Signal width differences
30 : /// - Amplification variations
31 : ///
32 : /// the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
33 : /// one chip and somewhat large between different chips.
34 : ///
35 : /// Histograms:
36 : /// For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum]) is created when
37 : /// it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
38 : /// TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
39 : ///
40 : ///
41 : /// Working principle:
42 : /// ------------------
43 : /// Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
44 : /// (see below). These in the end call the Update(...) function.
45 : ///
46 : /// - the Update(...) function:
47 : /// In this function the array fPadSignal is filled with the adc signals between the specified range
48 : /// fFirstTimeBin and fLastTimeBin for the current pad.
49 : /// before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
50 : /// stored in fPadSignal.
51 : ///
52 : /// - the ProcessPad() function:
53 : /// Find Pedestal and Noise information
54 : /// - use database information which has to be set by calling
55 : /// SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
56 : /// - if no information from the pedestal data base
57 : /// is available the informaion is calculated on the fly ( see FindPedestal() function )
58 : ///
59 : /// Find the Pulser signal information
60 : /// - calculate mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
61 : /// the Q sum is scaled by pad area
62 : /// (see FindPulserSignal(...) function)
63 : ///
64 : /// Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
65 : /// Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
66 : ///
67 : /// At the end of each event the EndEvent() function is called
68 : ///
69 : /// - the EndEvent() function:
70 : /// calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
71 : /// This is done to overcome syncronisation problems between the trigger and the fec clock.
72 : ///
73 : /// After accumulating the desired statistics the Analyse() function has to be called.
74 : /// - the Analyse() function
75 : /// Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
76 : /// the AliMathBase::GetCOG(...) function, and the calibration
77 : /// storage classes (AliTPCCalROC) are filled for each ROC.
78 : /// The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
79 : /// fCalRocArrayQ;
80 : ///
81 : ///
82 : ///
83 : /// User interface for filling data:
84 : /// --------------------------------
85 : ///
86 : /// To Fill information one of the following functions can be used:
87 : ///
88 : /// Bool_t ProcessEvent(eventHeaderStruct *event);
89 : /// - process Date event
90 : /// - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
91 : ///
92 : /// Bool_t ProcessEvent(AliRawReader *rawReader);
93 : /// - process AliRawReader event
94 : /// - use AliTPCRawStreamV3 to loop over data and call ProcessEvent(AliTPCRawStreamV3 *rawStream)
95 : ///
96 : /// Bool_t ProcessEvent(AliTPCRawStreamV3 *rawStream);
97 : /// - process event from AliTPCRawStreamV3
98 : /// - call Update function for signal filling
99 : ///
100 : /// Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
101 : /// iPad, const Int_t iTimeBin, const Float_t signal);
102 : /// - directly fill signal information (sector, row, pad, time bin, pad)
103 : /// to the reference histograms
104 : ///
105 : /// It is also possible to merge two independently taken calibrations using the function
106 : ///
107 : /// void Merge(AliTPCCalibPulser *sig)
108 : /// - copy histograms in 'sig' if the do not exist in this instance
109 : /// - Add histograms in 'sig' to the histograms in this instance if the allready exist
110 : /// - After merging call Analyse again!
111 : ///
112 : ///
113 : ///
114 : /// -- example: filling data using root raw data:
115 : /// ~~~{.cpp}
116 : /// void fillSignal(Char_t *filename)
117 : /// {
118 : /// rawReader = new AliRawReaderRoot(fileName);
119 : /// if ( !rawReader ) return;
120 : /// AliTPCCalibPulser *calib = new AliTPCCalibPulser;
121 : /// while (rawReader->NextEvent()){
122 : /// calib->ProcessEvent(rawReader);
123 : /// }
124 : /// calib->Analyse();
125 : /// calib->DumpToFile("SignalData.root");
126 : /// delete rawReader;
127 : /// delete calib;
128 : /// }
129 : /// ~~~
130 : ///
131 : /// What kind of information is stored and how to retrieve them:
132 : /// ------------------------------------------------------------
133 : ///
134 : /// - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
135 : ///
136 : /// TH2F *GetHistoT0(Int_t sector);
137 : /// TH2F *GetHistoRMS(Int_t sector);
138 : /// TH2F *GetHistoQ(Int_t sector);
139 : ///
140 : /// - Accessing the calibration storage objects:
141 : ///
142 : /// AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values
143 : /// AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values
144 : /// AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values
145 : ///
146 : /// example for visualisation:
147 : /// if the file "SignalData.root" was created using the above example one could do the following:
148 : ///
149 : /// ~~~{.cpp}
150 : /// TFile fileSignal("SignalData.root")
151 : /// AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
152 : /// sig->GetCalRocT0(0)->Draw("colz");
153 : /// sig->GetCalRocRMS(0)->Draw("colz");
154 : /// ~~~
155 : ///
156 : /// or use the AliTPCCalPad functionality:
157 : ///
158 : /// ~~~{.cpp}
159 : /// AliTPCCalPad padT0(ped->GetCalPadT0());
160 : /// AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
161 : /// padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
162 : /// padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
163 : /// ~~~
164 :
165 : //Root includes
166 : #include <TObjArray.h>
167 : #include <TH1F.h>
168 : #include <TH2F.h>
169 : #include <TH2S.h>
170 : #include <TString.h>
171 : #include <TVectorF.h>
172 : #include <TMath.h>
173 : #include <TMap.h>
174 : #include <TDirectory.h>
175 : #include <TSystem.h>
176 : #include <TROOT.h>
177 : #include <TFile.h>
178 :
179 : //AliRoot includes
180 : #include "AliRawReader.h"
181 : #include "AliRawReaderRoot.h"
182 : #include "AliRawReaderDate.h"
183 : #include "AliTPCCalROC.h"
184 : #include "AliTPCCalPad.h"
185 : #include "AliTPCROC.h"
186 : #include "AliTPCParam.h"
187 : #include "AliTPCCalibPulser.h"
188 : #include "AliTPCcalibDB.h"
189 : #include "AliMathBase.h"
190 : #include "AliLog.h"
191 : #include "TTreeStream.h"
192 :
193 : //date
194 : #include "event.h"
195 :
196 :
197 :
198 :
199 : /// \cond CLASSIMP
200 24 : ClassImp(AliTPCCalibPulser)
201 : /// \endcond
202 :
203 : AliTPCCalibPulser::AliTPCCalibPulser() :
204 0 : AliTPCCalibRawBase(),
205 0 : fNbinsT0(200),
206 0 : fXminT0(-2),
207 0 : fXmaxT0(2),
208 0 : fNbinsQ(200),
209 0 : fXminQ(10),
210 0 : fXmaxQ(40),
211 0 : fNbinsRMS(100),
212 0 : fXminRMS(0.1),
213 0 : fXmaxRMS(5.1),
214 0 : fPeakIntMinus(2),
215 0 : fPeakIntPlus(2),
216 0 : fIsZeroSuppressed(kFALSE),
217 0 : fLastSector(-1),
218 0 : fParam(new AliTPCParam),
219 0 : fPedestalTPC(0x0),
220 0 : fPadNoiseTPC(0x0),
221 0 : fOutliers(0x0),
222 0 : fPedestalROC(0x0),
223 0 : fPadNoiseROC(0x0),
224 0 : fCalRocArrayT0(72),
225 0 : fCalRocArrayQ(72),
226 0 : fCalRocArrayRMS(72),
227 0 : fCalRocArrayOutliers(72),
228 0 : fHistoQArray(72),
229 0 : fHistoT0Array(72),
230 0 : fHistoRMSArray(72),
231 0 : fHMeanTimeSector(0x0),
232 0 : fVMeanTimeSector(72),
233 0 : fPadTimesArrayEvent(72),
234 0 : fPadQArrayEvent(72),
235 0 : fPadRMSArrayEvent(72),
236 0 : fPadPedestalArrayEvent(72),
237 0 : fCurrentChannel(-1),
238 0 : fCurrentSector(-1),
239 0 : fCurrentRow(-1),
240 0 : fCurrentPad(-1),
241 0 : fMaxPadSignal(-1),
242 0 : fMaxTimeBin(-1),
243 0 : fPadSignal(1024),
244 0 : fPadPedestal(0),
245 0 : fPadNoise(0),
246 0 : fVTime0Offset(72),
247 0 : fVTime0OffsetCounter(72)
248 0 : {
249 : //
250 : // AliTPCSignal default constructor
251 : //
252 0 : SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
253 0 : fFirstTimeBin=60;
254 0 : fLastTimeBin=900;
255 0 : fParam->Update();
256 0 : }
257 : //_____________________________________________________________________
258 : AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
259 0 : AliTPCCalibRawBase(sig),
260 0 : fNbinsT0(sig.fNbinsT0),
261 0 : fXminT0(sig.fXminT0),
262 0 : fXmaxT0(sig.fXmaxT0),
263 0 : fNbinsQ(sig.fNbinsQ),
264 0 : fXminQ(sig.fXminQ),
265 0 : fXmaxQ(sig.fXmaxQ),
266 0 : fNbinsRMS(sig.fNbinsRMS),
267 0 : fXminRMS(sig.fXminRMS),
268 0 : fXmaxRMS(sig.fXmaxRMS),
269 0 : fPeakIntMinus(sig.fPeakIntMinus),
270 0 : fPeakIntPlus(sig.fPeakIntPlus),
271 0 : fIsZeroSuppressed(sig.fIsZeroSuppressed),
272 0 : fLastSector(-1),
273 0 : fParam(new AliTPCParam),
274 0 : fPedestalTPC(0x0),
275 0 : fPadNoiseTPC(0x0),
276 0 : fOutliers(0x0),
277 0 : fPedestalROC(0x0),
278 0 : fPadNoiseROC(0x0),
279 0 : fCalRocArrayT0(72),
280 0 : fCalRocArrayQ(72),
281 0 : fCalRocArrayRMS(72),
282 0 : fCalRocArrayOutliers(72),
283 0 : fHistoQArray(72),
284 0 : fHistoT0Array(72),
285 0 : fHistoRMSArray(72),
286 0 : fHMeanTimeSector(0x0),
287 0 : fVMeanTimeSector(72),
288 0 : fPadTimesArrayEvent(72),
289 0 : fPadQArrayEvent(72),
290 0 : fPadRMSArrayEvent(72),
291 0 : fPadPedestalArrayEvent(72),
292 0 : fCurrentChannel(-1),
293 0 : fCurrentSector(-1),
294 0 : fCurrentRow(-1),
295 0 : fCurrentPad(-1),
296 0 : fMaxPadSignal(-1),
297 0 : fMaxTimeBin(-1),
298 0 : fPadSignal(1024),
299 0 : fPadPedestal(0),
300 0 : fPadNoise(0),
301 0 : fVTime0Offset(72),
302 0 : fVTime0OffsetCounter(72)
303 0 : {
304 : /// AliTPCSignal default constructor
305 :
306 0 : for (Int_t iSec = 0; iSec < 72; ++iSec){
307 0 : const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
308 0 : const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
309 0 : const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
310 0 : const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
311 :
312 0 : const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
313 0 : const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
314 0 : const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
315 :
316 0 : if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
317 0 : if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
318 0 : if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
319 0 : if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
320 :
321 0 : if ( hQ != 0x0 ){
322 0 : TH2S *hNew = new TH2S(*hQ);
323 0 : hNew->SetDirectory(0);
324 0 : fHistoQArray.AddAt(hNew,iSec);
325 0 : }
326 0 : if ( hT0 != 0x0 ){
327 0 : TH2S *hNew = new TH2S(*hT0);
328 0 : hNew->SetDirectory(0);
329 0 : fHistoQArray.AddAt(hNew,iSec);
330 0 : }
331 0 : if ( hRMS != 0x0 ){
332 0 : TH2S *hNew = new TH2S(*hRMS);
333 0 : hNew->SetDirectory(0);
334 0 : fHistoQArray.AddAt(hNew,iSec);
335 0 : }
336 0 : fVMeanTimeSector[iSec]=sig.fVMeanTimeSector[iSec];
337 : }
338 :
339 0 : if ( sig.fHMeanTimeSector ) fHMeanTimeSector=(TH2F*)sig.fHMeanTimeSector->Clone();
340 0 : fParam->Update();
341 0 : }
342 : AliTPCCalibPulser::AliTPCCalibPulser(const TMap *config) :
343 0 : AliTPCCalibRawBase(),
344 0 : fNbinsT0(200),
345 0 : fXminT0(-2),
346 0 : fXmaxT0(2),
347 0 : fNbinsQ(200),
348 0 : fXminQ(10),
349 0 : fXmaxQ(40),
350 0 : fNbinsRMS(100),
351 0 : fXminRMS(0.1),
352 0 : fXmaxRMS(5.1),
353 0 : fPeakIntMinus(2),
354 0 : fPeakIntPlus(2),
355 0 : fIsZeroSuppressed(kFALSE),
356 0 : fLastSector(-1),
357 0 : fParam(new AliTPCParam),
358 0 : fPedestalTPC(0x0),
359 0 : fPadNoiseTPC(0x0),
360 0 : fOutliers(0x0),
361 0 : fPedestalROC(0x0),
362 0 : fPadNoiseROC(0x0),
363 0 : fCalRocArrayT0(72),
364 0 : fCalRocArrayQ(72),
365 0 : fCalRocArrayRMS(72),
366 0 : fCalRocArrayOutliers(72),
367 0 : fHistoQArray(72),
368 0 : fHistoT0Array(72),
369 0 : fHistoRMSArray(72),
370 0 : fHMeanTimeSector(0x0),
371 0 : fVMeanTimeSector(72),
372 0 : fPadTimesArrayEvent(72),
373 0 : fPadQArrayEvent(72),
374 0 : fPadRMSArrayEvent(72),
375 0 : fPadPedestalArrayEvent(72),
376 0 : fCurrentChannel(-1),
377 0 : fCurrentSector(-1),
378 0 : fCurrentRow(-1),
379 0 : fCurrentPad(-1),
380 0 : fMaxPadSignal(-1),
381 0 : fMaxTimeBin(-1),
382 0 : fPadSignal(1024),
383 0 : fPadPedestal(0),
384 0 : fPadNoise(0),
385 0 : fVTime0Offset(72),
386 0 : fVTime0OffsetCounter(72)
387 0 : {
388 : /// This constructor uses a TMap for setting some parametes
389 :
390 0 : SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
391 0 : fFirstTimeBin=60;
392 0 : fLastTimeBin=900;
393 0 : if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
394 0 : if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
395 0 : if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
396 0 : if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
397 0 : if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
398 0 : if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
399 0 : if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
400 0 : if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
401 0 : if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
402 0 : if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
403 0 : if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
404 0 : if (config->GetValue("PeakIntMinus")) fPeakIntMinus = (Int_t)((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atof();
405 0 : if (config->GetValue("PeakIntPlus")) fPeakIntPlus = (Int_t)((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atof();
406 0 : if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
407 :
408 0 : fParam->Update();
409 0 : }
410 : //_____________________________________________________________________
411 : AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
412 : {
413 : /// assignment operator
414 :
415 0 : if (&source == this) return *this;
416 0 : new (this) AliTPCCalibPulser(source);
417 :
418 0 : return *this;
419 0 : }
420 : //_____________________________________________________________________
421 : AliTPCCalibPulser::~AliTPCCalibPulser()
422 0 : {
423 : /// destructor
424 :
425 0 : Reset();
426 0 : delete fParam;
427 0 : }
428 : void AliTPCCalibPulser::Reset()
429 : {
430 : /// Delete all information: Arrays, Histograms, CalRoc objects
431 :
432 0 : fCalRocArrayT0.Delete();
433 0 : fCalRocArrayQ.Delete();
434 0 : fCalRocArrayRMS.Delete();
435 0 : fCalRocArrayOutliers.Delete();
436 :
437 0 : fHistoQArray.Delete();
438 0 : fHistoT0Array.Delete();
439 0 : fHistoRMSArray.Delete();
440 :
441 0 : fPadTimesArrayEvent.Delete();
442 0 : fPadQArrayEvent.Delete();
443 0 : fPadRMSArrayEvent.Delete();
444 0 : fPadPedestalArrayEvent.Delete();
445 :
446 0 : if (fHMeanTimeSector) delete fHMeanTimeSector;
447 0 : }
448 : //_____________________________________________________________________
449 : Int_t AliTPCCalibPulser::Update(const Int_t icsector,
450 : const Int_t icRow,
451 : const Int_t icPad,
452 : const Int_t icTimeBin,
453 : const Float_t csignal)
454 : {
455 : /// Signal filling methode on the fly pedestal and time offset correction if necessary.
456 : /// no extra analysis necessary. Assumes knowledge of the signal shape!
457 : /// assumes that it is looped over consecutive time bins of one pad
458 :
459 0 : if (icRow<0) return 0;
460 0 : if (icPad<0) return 0;
461 0 : if (icTimeBin<0) return 0;
462 0 : if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
463 :
464 0 : if ( icRow<0 || icPad<0 ){
465 0 : AliWarning("Wrong Pad or Row number, skipping!");
466 0 : return 0;
467 : }
468 :
469 0 : Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
470 :
471 : //init first pad and sector in this event
472 0 : if ( fCurrentChannel == -1 ) {
473 0 : fCurrentChannel = iChannel;
474 0 : fCurrentSector = icsector;
475 0 : fCurrentRow = icRow;
476 0 : fCurrentPad = icPad;
477 0 : }
478 :
479 : //process last pad if we change to a new one
480 0 : if ( iChannel != fCurrentChannel ){
481 0 : ProcessPad();
482 0 : fLastSector=fCurrentSector;
483 0 : fCurrentChannel = iChannel;
484 0 : fCurrentSector = icsector;
485 0 : fCurrentRow = icRow;
486 0 : fCurrentPad = icPad;
487 0 : }
488 :
489 : //fill signals for current pad
490 0 : fPadSignal[icTimeBin]=csignal;
491 0 : if ( csignal > fMaxPadSignal ){
492 0 : fMaxPadSignal = csignal;
493 0 : fMaxTimeBin = icTimeBin;
494 0 : }
495 : return 0;
496 0 : }
497 : //_____________________________________________________________________
498 : void AliTPCCalibPulser::FindPedestal(Float_t part)
499 : {
500 : /// find pedestal and noise for the current pad. Use either database or
501 : /// truncated mean with part*100%
502 :
503 : Bool_t noPedestal = kTRUE;;
504 0 : if (fPedestalTPC&&fPadNoiseTPC){
505 : //use pedestal database
506 : //only load new pedestals if the sector has changed
507 0 : if ( fCurrentSector!=fLastSector ){
508 0 : fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
509 0 : fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
510 0 : }
511 :
512 0 : if ( fPedestalROC&&fPadNoiseROC ){
513 0 : fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
514 0 : fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
515 : noPedestal = kFALSE;
516 0 : }
517 :
518 : }
519 :
520 : //if we are not running with pedestal database, or for the current sector there is no information
521 : //available, calculate the pedestal and noise on the fly
522 0 : if ( noPedestal ) {
523 : const Int_t kPedMax = 100; //maximum pedestal value
524 : Float_t max = 0;
525 : Int_t median = -1;
526 : Int_t count0 = 0;
527 : Int_t count1 = 0;
528 : //
529 : Float_t padSignal=0;
530 : //
531 0 : UShort_t histo[kPedMax];
532 0 : memset(histo,0,kPedMax*sizeof(UShort_t));
533 :
534 0 : for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
535 0 : padSignal = fPadSignal.GetMatrixArray()[i];
536 0 : if (padSignal<=0) continue;
537 0 : if (padSignal>max && i>10) {
538 : max = padSignal;
539 0 : }
540 0 : if (padSignal>kPedMax-1) continue;
541 0 : histo[Int_t(padSignal+0.5)]++;
542 0 : count0++;
543 0 : }
544 : //
545 0 : for (Int_t i=1; i<kPedMax; ++i){
546 0 : if (count1<count0*0.5) median=i;
547 0 : count1+=histo[i];
548 : }
549 : // truncated mean
550 : //
551 : // what if by chance histo[median] == 0 ?!?
552 0 : Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
553 : //
554 0 : for (Int_t idelta=1; idelta<10; ++idelta){
555 0 : if (median-idelta<=0) continue;
556 0 : if (median+idelta>kPedMax) continue;
557 0 : if (count<part*count1){
558 0 : count+=histo[median-idelta];
559 0 : mean +=histo[median-idelta]*(median-idelta);
560 0 : rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
561 0 : count+=histo[median+idelta];
562 0 : mean +=histo[median+idelta]*(median+idelta);
563 0 : rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
564 0 : }
565 : }
566 0 : fPadPedestal = 0;
567 0 : fPadNoise = 0;
568 0 : if ( count > 0 ) {
569 0 : mean/=count;
570 0 : rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
571 0 : fPadPedestal = mean;
572 0 : fPadNoise = rms;
573 0 : }
574 0 : }
575 0 : fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
576 0 : }
577 : //_____________________________________________________________________
578 : void AliTPCCalibPulser::FindPulserSignal(TVectorD ¶m, Float_t &qSum)
579 : {
580 : /// Find position, signal width and height of the CE signal (last signal)
581 : /// param[0] = Qmax, param[1] = mean time, param[2] = rms;
582 : /// maxima: array of local maxima of the pad signal use the one closest to the mean CE position
583 :
584 : Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
585 0 : Int_t cemaxpos = fMaxTimeBin;
586 0 : Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal sum
587 0 : Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal max
588 0 : const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
589 0 : const Int_t kCemax = fPeakIntPlus;
590 0 : param[0] = ceQmax;
591 0 : param[1] = ceTime;
592 0 : param[2] = ceRMS;
593 0 : qSum = ceQsum;
594 :
595 0 : if (cemaxpos>0){
596 0 : ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
597 0 : if ( ceQmax<ceMaxThreshold ) return;
598 0 : for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
599 0 : Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
600 0 : if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
601 0 : ceTime+=signal*(i+0.5);
602 0 : ceRMS +=signal*(i+0.5)*(i+0.5);
603 0 : ceQsum+=signal;
604 0 : }
605 : }
606 0 : }
607 0 : if (ceQsum>ceSumThreshold) {
608 0 : ceTime/=ceQsum;
609 0 : ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
610 0 : ceTime-=GetL1PhaseTB();
611 : //only fill the Time0Offset if pad was not marked as an outlier!
612 0 : if ( !fOutliers ){
613 : //skip edge pads for calculating the mean time
614 0 : if ( !IsEdgePad(fCurrentSector,fCurrentRow,fCurrentPad) ){
615 0 : fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
616 0 : fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
617 0 : GetHistoTSec()->Fill(ceTime,fCurrentSector);
618 0 : }
619 : } else {
620 0 : if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
621 0 : fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
622 0 : fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
623 0 : }
624 : }
625 :
626 : //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
627 : // the pick-up signal should scale with the pad area. In addition
628 : // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
629 : // ratio 2/3. The pad area we express in mm2 (factor 100). We normalise the signal
630 : // to the OROC signal (factor 2/3 for the IROCs).
631 0 : Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow)*100;
632 0 : if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
633 :
634 0 : ceQsum/=norm;
635 0 : } else {
636 : ceQmax=0;
637 : ceTime=0;
638 : ceRMS =0;
639 : ceQsum=0;
640 : }
641 0 : param[0] = ceQmax;
642 0 : param[1] = ceTime;
643 0 : param[2] = ceRMS;
644 0 : qSum = ceQsum;
645 0 : }
646 : //_____________________________________________________________________
647 : void AliTPCCalibPulser::ProcessPad()
648 : {
649 : /// Process data of current pad
650 :
651 0 : FindPedestal();
652 0 : TVectorD param(3);
653 0 : Float_t qSum;
654 0 : FindPulserSignal(param, qSum);
655 :
656 0 : Double_t meanT = param[1];
657 0 : Double_t sigmaT = param[2];
658 :
659 :
660 : //Fill Event T0 counter
661 0 : (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
662 :
663 : //Fill Q histogram
664 : // GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel ); //use linear scale, needs also a change in the Analyse funciton.
665 0 : GetHistoQ(fCurrentSector,kTRUE)->Fill( qSum, fCurrentChannel );
666 :
667 : //Fill RMS histogram
668 0 : GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
669 :
670 :
671 : //Fill debugging info
672 0 : if ( GetStreamLevel()>0 ){
673 0 : TTreeSRedirector *streamer=GetDebugStreamer();
674 0 : if ( GetStreamLevel() == 1 ){
675 0 : if ( streamer ) {
676 0 : Int_t padc = fCurrentPad-(fROC->GetNPads(fCurrentSector,fCurrentRow)/2);
677 0 : (*streamer) << "PadSignals" <<
678 0 : "Event=" <<fNevents <<
679 0 : "Sector=" <<fCurrentSector<<
680 0 : "Row=" <<fCurrentRow<<
681 0 : "Pad=" <<fCurrentPad<<
682 0 : "PadC=" <<padc<<
683 0 : "Channel="<<fCurrentChannel<<
684 0 : "Sum=" <<qSum<<
685 0 : "params.="<<¶m<<
686 0 : "signal.=" <<&fPadSignal<<
687 : "\n";
688 0 : }
689 : } else { //debug > 1
690 0 : (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
691 0 : (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
692 0 : (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
693 : }
694 0 : }
695 0 : ResetPad();
696 0 : }
697 : //_____________________________________________________________________
698 : void AliTPCCalibPulser::EndEvent()
699 : {
700 : /// Process data of current event
701 :
702 : //check if last pad has allready been processed, if not do so
703 0 : if ( fMaxTimeBin>-1 ) ProcessPad();
704 :
705 : //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
706 0 : for ( Int_t iSec = 0; iSec<72; ++iSec ){
707 0 : TVectorF *vTimes = GetPadTimesEvent(iSec);
708 0 : if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
709 0 : Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
710 0 : for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
711 0 : Float_t time = (*vTimes).GetMatrixArray()[iChannel];
712 :
713 0 : GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
714 : //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
715 :
716 :
717 : //Debug start
718 0 : if ( GetStreamLevel()>1 ){
719 0 : TTreeSRedirector *streamer=GetDebugStreamer();
720 0 : if ( streamer ) {
721 0 : Int_t row=0;
722 0 : Int_t pad=0;
723 0 : Int_t padc=0;
724 :
725 0 : Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
726 0 : Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
727 :
728 0 : UInt_t channel=iChannel;
729 0 : Int_t sector=iSec;
730 :
731 0 : while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
732 0 : pad = channel-fROC->GetRowIndexes(sector)[row];
733 0 : padc = pad-(fROC->GetNPads(sector,row)/2);
734 :
735 0 : (*streamer) << "DataPad" <<
736 0 : "Event=" << fNevents <<
737 0 : "Sector="<< sector <<
738 0 : "Row=" << row<<
739 0 : "Pad=" << pad <<
740 0 : "PadC=" << padc <<
741 0 : "PadSec="<< channel <<
742 0 : "Time0=" << time0 <<
743 0 : "Time=" << time <<
744 0 : "RMS=" << rms <<
745 0 : "Sum=" << q <<
746 : "\n";
747 0 : }
748 0 : }
749 : //Debug end
750 0 : }
751 0 : }
752 0 : ++fNevents;
753 0 : }
754 : //_____________________________________________________________________
755 : TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
756 : Int_t nbinsY, Float_t ymin, Float_t ymax,
757 : const Char_t *type, Bool_t force)
758 : {
759 : /// return pointer to Q histogram
760 : /// if force is true create a new histogram if it doesn't exist allready
761 :
762 0 : if ( !force || arr->UncheckedAt(sector) )
763 0 : return (TH2S*)arr->UncheckedAt(sector);
764 :
765 : // if we are forced and histogram doesn't yes exist create it
766 : // new histogram with Q calib information. One value for each pad!
767 0 : TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
768 0 : nbinsY, ymin, ymax,
769 0 : fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
770 0 : hist->SetDirectory(0);
771 0 : arr->AddAt(hist,sector);
772 : return hist;
773 0 : }
774 : //_____________________________________________________________________
775 : TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
776 : {
777 : /// return pointer to T0 histogram
778 : /// if force is true create a new histogram if it doesn't exist allready
779 :
780 0 : TObjArray *arr = &fHistoT0Array;
781 0 : return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
782 : }
783 : //_____________________________________________________________________
784 : TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
785 : {
786 : /// return pointer to Q histogram
787 : /// if force is true create a new histogram if it doesn't exist allready
788 :
789 0 : TObjArray *arr = &fHistoQArray;
790 0 : return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
791 : }
792 : //_____________________________________________________________________
793 : TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
794 : {
795 : /// return pointer to Q histogram
796 : /// if force is true create a new histogram if it doesn't exist allready
797 :
798 0 : TObjArray *arr = &fHistoRMSArray;
799 0 : return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
800 : }
801 : //_____________________________________________________________________
802 : TH2F* AliTPCCalibPulser::GetHistoTSec()
803 : {
804 : /// return the pointer to the abs time distribution per sector
805 : /// create it if it does not exist
806 :
807 0 : if ( !fHMeanTimeSector ) //!!!if you change the binning here, you should also change it in the Analyse Function!!
808 0 : fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
809 0 : 20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
810 : 72,0,72);
811 0 : return fHMeanTimeSector;
812 0 : }
813 : //_____________________________________________________________________
814 : TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
815 : {
816 : /// return pointer to Pad Info from 'arr' for the current event and sector
817 : /// if force is true create it if it doesn't exist allready
818 :
819 0 : if ( !force || arr->UncheckedAt(sector) )
820 0 : return (TVectorF*)arr->UncheckedAt(sector);
821 :
822 0 : TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
823 0 : arr->AddAt(vect,sector);
824 : return vect;
825 0 : }
826 : //_____________________________________________________________________
827 : TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
828 : {
829 : /// return pointer to Pad Times Array for the current event and sector
830 : /// if force is true create it if it doesn't exist allready
831 :
832 0 : TObjArray *arr = &fPadTimesArrayEvent;
833 0 : return GetPadInfoEvent(sector,arr,force);
834 : }
835 : //_____________________________________________________________________
836 : TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
837 : {
838 : /// return pointer to Pad Q Array for the current event and sector
839 : /// if force is true create it if it doesn't exist allready
840 : /// for debugging purposes only
841 :
842 0 : TObjArray *arr = &fPadQArrayEvent;
843 0 : return GetPadInfoEvent(sector,arr,force);
844 : }
845 : //_____________________________________________________________________
846 : TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
847 : {
848 : /// return pointer to Pad RMS Array for the current event and sector
849 : /// if force is true create it if it doesn't exist allready
850 : /// for debugging purposes only
851 :
852 0 : TObjArray *arr = &fPadRMSArrayEvent;
853 0 : return GetPadInfoEvent(sector,arr,force);
854 : }
855 : //_____________________________________________________________________
856 : TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
857 : {
858 : /// return pointer to Pad RMS Array for the current event and sector
859 : /// if force is true create it if it doesn't exist allready
860 : /// for debugging purposes only
861 :
862 0 : TObjArray *arr = &fPadPedestalArrayEvent;
863 0 : return GetPadInfoEvent(sector,arr,force);
864 : }
865 : //_____________________________________________________________________
866 : AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
867 : {
868 : /// return pointer to ROC Calibration
869 : /// if force is true create a new histogram if it doesn't exist allready
870 :
871 0 : if ( !force || arr->UncheckedAt(sector) )
872 0 : return (AliTPCCalROC*)arr->UncheckedAt(sector);
873 :
874 : // if we are forced and histogram doesn't yes exist create it
875 :
876 : // new AliTPCCalROC for T0 information. One value for each pad!
877 0 : AliTPCCalROC *croc = new AliTPCCalROC(sector);
878 0 : arr->AddAt(croc,sector);
879 : return croc;
880 0 : }
881 : //_____________________________________________________________________
882 : AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
883 : {
884 : /// return pointer to Carge ROC Calibration
885 : /// if force is true create a new histogram if it doesn't exist allready
886 :
887 0 : TObjArray *arr = &fCalRocArrayT0;
888 0 : return GetCalRoc(sector, arr, force);
889 : }
890 : //_____________________________________________________________________
891 : AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
892 : {
893 : /// return pointer to T0 ROC Calibration
894 : /// if force is true create a new histogram if it doesn't exist allready
895 :
896 0 : TObjArray *arr = &fCalRocArrayQ;
897 0 : return GetCalRoc(sector, arr, force);
898 : }
899 : //_____________________________________________________________________
900 : AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
901 : {
902 : /// return pointer to signal width ROC Calibration
903 : /// if force is true create a new histogram if it doesn't exist allready
904 :
905 0 : TObjArray *arr = &fCalRocArrayRMS;
906 0 : return GetCalRoc(sector, arr, force);
907 : }
908 : //_____________________________________________________________________
909 : AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
910 : {
911 : /// return pointer to Outliers
912 : /// if force is true create a new histogram if it doesn't exist allready
913 :
914 0 : TObjArray *arr = &fCalRocArrayOutliers;
915 0 : return GetCalRoc(sector, arr, force);
916 : }
917 : //_____________________________________________________________________
918 : void AliTPCCalibPulser::ResetEvent()
919 : {
920 : /// Reset global counters -- Should be called before each event is processed
921 :
922 0 : fLastSector=-1;
923 0 : fCurrentSector=-1;
924 0 : fCurrentRow=-1;
925 0 : fCurrentPad=-1;
926 0 : fCurrentChannel=-1;
927 :
928 0 : ResetPad();
929 :
930 0 : fPadTimesArrayEvent.Delete();
931 :
932 0 : fPadQArrayEvent.Delete();
933 0 : fPadRMSArrayEvent.Delete();
934 0 : fPadPedestalArrayEvent.Delete();
935 :
936 0 : for ( Int_t i=0; i<72; ++i ){
937 0 : fVTime0Offset[i]=0;
938 0 : fVTime0OffsetCounter[i]=0;
939 : }
940 0 : }
941 : //_____________________________________________________________________
942 : void AliTPCCalibPulser::ResetPad()
943 : {
944 : /// Reset pad infos -- Should be called after a pad has been processed
945 :
946 0 : for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
947 0 : fPadSignal[i] = 0;
948 0 : fMaxTimeBin = -1;
949 0 : fMaxPadSignal = -1;
950 0 : fPadPedestal = -1;
951 0 : fPadNoise = -1;
952 0 : }
953 : //_____________________________________________________________________
954 : Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
955 : {
956 : /// return true if pad is on the edge of a row
957 :
958 : Int_t edge1 = 0;
959 0 : Int_t edge2 = fROC->GetNPads(sector,row)-1;
960 0 : if ( pad == edge1 || pad == edge2 ) return kTRUE;
961 :
962 0 : return kFALSE;
963 0 : }
964 : //_____________________________________________________________________
965 : void AliTPCCalibPulser::Merge(AliTPCCalibPulser * const sig)
966 : {
967 : /// Merge reference histograms of sig to the current AliTPCCalibPulser
968 :
969 0 : MergeBase(sig);
970 : //merge histograms
971 0 : for (Int_t iSec=0; iSec<72; ++iSec){
972 0 : TH2S *hRefQmerge = sig->GetHistoQ(iSec);
973 0 : TH2S *hRefT0merge = sig->GetHistoT0(iSec);
974 0 : TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
975 :
976 :
977 0 : if ( hRefQmerge ){
978 0 : TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
979 0 : TH2S *hRefQ = GetHistoQ(iSec);
980 0 : if ( hRefQ ) hRefQ->Add(hRefQmerge);
981 : else {
982 0 : TH2S *hist = new TH2S(*hRefQmerge);
983 0 : hist->SetDirectory(0);
984 0 : fHistoQArray.AddAt(hist, iSec);
985 : }
986 0 : hRefQmerge->SetDirectory(dir);
987 0 : }
988 0 : if ( hRefT0merge ){
989 0 : TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
990 0 : TH2S *hRefT0 = GetHistoT0(iSec);
991 0 : if ( hRefT0 ) hRefT0->Add(hRefT0merge);
992 : else {
993 0 : TH2S *hist = new TH2S(*hRefT0merge);
994 0 : hist->SetDirectory(0);
995 0 : fHistoT0Array.AddAt(hist, iSec);
996 : }
997 0 : hRefT0merge->SetDirectory(dir);
998 0 : }
999 0 : if ( hRefRMSmerge ){
1000 0 : TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1001 0 : TH2S *hRefRMS = GetHistoRMS(iSec);
1002 0 : if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1003 : else {
1004 0 : TH2S *hist = new TH2S(*hRefRMSmerge);
1005 0 : hist->SetDirectory(0);
1006 0 : fHistoRMSArray.AddAt(hist, iSec);
1007 : }
1008 0 : hRefRMSmerge->SetDirectory(dir);
1009 0 : }
1010 :
1011 : }
1012 0 : if ( sig->fHMeanTimeSector ){
1013 0 : TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1014 0 : if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1015 : else {
1016 0 : fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1017 0 : fHMeanTimeSector->SetDirectory(0);
1018 : }
1019 0 : sig->fHMeanTimeSector->SetDirectory(dir);
1020 0 : }
1021 0 : }
1022 :
1023 :
1024 : //_____________________________________________________________________
1025 : Long64_t AliTPCCalibPulser::Merge(TCollection * const list)
1026 : {
1027 : /// Merge all objects of this type in list
1028 :
1029 : Long64_t nmerged=1;
1030 :
1031 0 : TIter next(list);
1032 : AliTPCCalibPulser *ce=0;
1033 : TObject *o=0;
1034 :
1035 0 : while ( (o=next()) ){
1036 0 : ce=dynamic_cast<AliTPCCalibPulser*>(o);
1037 0 : if (ce){
1038 0 : Merge(ce);
1039 0 : ++nmerged;
1040 0 : }
1041 : }
1042 :
1043 : return nmerged;
1044 0 : }
1045 :
1046 : //_____________________________________________________________________
1047 : void AliTPCCalibPulser::Analyse()
1048 : {
1049 : /// Calculate calibration constants
1050 :
1051 0 : TVectorD paramQ(3);
1052 0 : TVectorD paramT0(3);
1053 0 : TVectorD paramRMS(3);
1054 0 : TMatrixD dummy(3,3);
1055 : //calculate mean time for each sector and mean time for each side
1056 0 : TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1057 0 : fVMeanTimeSector.Zero();
1058 :
1059 0 : for (Int_t iSec=0; iSec<72; ++iSec){
1060 0 : TH2S *hT0 = GetHistoT0(iSec);
1061 0 : if (!hT0 ) continue;
1062 : //calculate sector mean T
1063 0 : if ( fHMeanTimeSector ){
1064 0 : Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1065 0 : Int_t offset = (nbinsT+2)*(iSec+1);
1066 0 : Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1067 : Int_t entries=0;
1068 0 : for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1069 0 : hMeanTsec.Set(nbinsT+2,arrP);
1070 0 : hMeanTsec.SetEntries(entries);
1071 0 : paramT0.Zero();
1072 : // truncated mean: remove lower 5% and upper 5%
1073 0 : if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,¶mT0,0.05,.95);
1074 0 : fVMeanTimeSector[iSec]=paramT0[1];
1075 0 : }
1076 :
1077 0 : AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1078 0 : AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1079 0 : AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1080 0 : AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1081 :
1082 0 : TH2S *hQ = GetHistoQ(iSec);
1083 0 : TH2S *hRMS = GetHistoRMS(iSec);
1084 :
1085 0 : Short_t *arrayhQ = hQ->GetArray();
1086 0 : Short_t *arrayhT0 = hT0->GetArray();
1087 0 : Short_t *arrayhRMS = hRMS->GetArray();
1088 :
1089 0 : UInt_t nChannels = fROC->GetNChannels(iSec);
1090 0 : Float_t meanTsec = fVMeanTimeSector[iSec];
1091 :
1092 : //debug
1093 0 : Int_t row=0;
1094 0 : Int_t pad=0;
1095 0 : Int_t padc=0;
1096 : //! debug
1097 :
1098 0 : for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1099 :
1100 0 : Float_t cogTime0 = -1000;
1101 0 : Float_t cogQ = -1000;
1102 0 : Float_t cogRMS = -1000;
1103 : Float_t cogOut = 0;
1104 :
1105 0 : Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1106 0 : Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1107 0 : Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1108 : /*
1109 : AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,¶mQ,&dummy);
1110 : AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,¶mT0,&dummy);
1111 : AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,¶mRMS,&dummy);
1112 : cogQ = paramQ[1];
1113 : cogTime0 = paramT0[1];
1114 : cogRMS = paramRMS[1];
1115 : */
1116 0 : cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1117 0 : cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1118 0 : cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1119 :
1120 : /*
1121 : if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1122 : cogOut = 1;
1123 : cogTime0 = 0;
1124 : cogQ = 0;
1125 : cogRMS = 0;
1126 : }
1127 : */
1128 : // rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1129 0 : rocQ->SetValue(iChannel, cogQ);
1130 0 : rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1131 0 : rocRMS->SetValue(iChannel, cogRMS);
1132 0 : rocOut->SetValue(iChannel, cogOut);
1133 :
1134 : // in case a channel has no data set the value to 0
1135 0 : if (TMath::Abs(cogTime0-fXminT0)<1e-10){
1136 0 : rocQ->SetValue(iChannel, 0);
1137 0 : rocT0->SetValue(iChannel, 0); //offset by mean time of the sector
1138 0 : rocRMS->SetValue(iChannel, 0);
1139 0 : }
1140 :
1141 : //debug
1142 0 : if ( GetStreamLevel() > 2 ){
1143 0 : TTreeSRedirector *streamer=GetDebugStreamer();
1144 0 : if ( streamer ) {
1145 0 : while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1146 0 : pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1147 0 : padc = pad-(fROC->GetNPads(iSec,row)/2);
1148 :
1149 0 : (*streamer) << "DataEnd" <<
1150 0 : "Sector=" << iSec <<
1151 0 : "Pad=" << pad <<
1152 0 : "PadC=" << padc <<
1153 0 : "Row=" << row <<
1154 0 : "PadSec=" << iChannel <<
1155 0 : "Q=" << cogQ <<
1156 0 : "T0=" << cogTime0 <<
1157 0 : "RMS=" << cogRMS <<
1158 : "\n";
1159 : }
1160 0 : }
1161 : //! debug
1162 0 : }
1163 :
1164 :
1165 0 : }
1166 0 : }
1167 : //_____________________________________________________________________
1168 : //_________________________ Test Functions ___________________________
1169 : //_____________________________________________________________________
1170 : TObjArray* AliTPCCalibPulser::TestBinning()
1171 : {
1172 : /// Function to test the binning of the reference histograms
1173 : /// type: T0, Q or RMS
1174 : /// mode: 0 - number of filled bins per channel
1175 : /// 1 - number of empty bins between filled bins in one ROC
1176 : /// returns TObjArray with the test histograms type*2+mode:
1177 : /// position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1178 :
1179 :
1180 0 : TObjArray *histArray = new TObjArray(6);
1181 0 : const Char_t *type[] = {"T0","Q","RMS"};
1182 0 : Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1183 :
1184 0 : for (Int_t itype = 0; itype<3; ++itype){
1185 0 : for (Int_t imode=0; imode<2; ++imode){
1186 0 : Int_t icount = itype*2+imode;
1187 0 : histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1188 0 : Form("Test Binning of '%s', mode - %d",type[itype],imode),
1189 : 72,0,72),
1190 : icount);
1191 : }
1192 : }
1193 :
1194 :
1195 : TH2S *hRef=0x0;
1196 : Short_t *array=0x0;
1197 0 : for (Int_t itype = 0; itype<3; ++itype){
1198 0 : for (Int_t iSec=0; iSec<72; ++iSec){
1199 0 : if ( itype == 0 ) hRef = GetHistoT0(iSec);
1200 0 : if ( itype == 1 ) hRef = GetHistoQ(iSec);
1201 0 : if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1202 0 : if ( hRef == 0x0 ) continue;
1203 0 : array = (hRef->GetArray());
1204 0 : UInt_t nChannels = fROC->GetNChannels(iSec);
1205 :
1206 : Int_t nempty=0;
1207 0 : for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1208 : Int_t nfilled=0;
1209 0 : Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1210 : Int_t c1 = 0;
1211 : Int_t c2 = 0;
1212 0 : for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1213 0 : if ( array[offset+iBin]>0 ) {
1214 0 : nfilled++;
1215 0 : if ( c1 && c2 ) nempty++;
1216 : else c1 = 1;
1217 : }
1218 0 : else if ( c1 ) c2 = 1;
1219 :
1220 : }
1221 0 : ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1222 : }
1223 0 : ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1224 0 : }
1225 : }
1226 0 : return histArray;
1227 0 : }
|
