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