]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Base/AliTPCCalibPulser.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
o remove printf
[u/mrichter/AliRoot.git] / TPC / Base / 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
829455ad 99 - use AliTPCRawStreamV3 to loop over data and call ProcessEvent(AliTPCRawStreamV3 *rawStream)
467a4b25 100
829455ad 101 Bool_t ProcessEvent(AliTPCRawStreamV3 *rawStream);
102 - process event from AliTPCRawStreamV3
467a4b25 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"
3cd27a08 184#include "AliTPCCalROC.h"
185#include "AliTPCCalPad.h"
186#include "AliTPCROC.h"
187#include "AliTPCParam.h"
188#include "AliTPCCalibPulser.h"
189#include "AliTPCcalibDB.h"
190#include "AliMathBase.h"
d8faeea7 191#include "AliLog.h"
3cd27a08 192#include "TTreeStream.h"
467a4b25 193
3cd27a08 194//date
195#include "event.h"
196
197
198
199
200ClassImp(AliTPCCalibPulser)
201
202AliTPCCalibPulser::AliTPCCalibPulser() :
c3066940 203AliTPCCalibRawBase(),
204fNbinsT0(200),
205fXminT0(-2),
206fXmaxT0(2),
207fNbinsQ(200),
208fXminQ(10),
209fXmaxQ(40),
210fNbinsRMS(100),
211fXminRMS(0.1),
212fXmaxRMS(5.1),
213fPeakIntMinus(2),
214fPeakIntPlus(2),
215fIsZeroSuppressed(kFALSE),
216fLastSector(-1),
217fParam(new AliTPCParam),
218fPedestalTPC(0x0),
219fPadNoiseTPC(0x0),
220fOutliers(0x0),
221fPedestalROC(0x0),
222fPadNoiseROC(0x0),
223fCalRocArrayT0(72),
224fCalRocArrayQ(72),
225fCalRocArrayRMS(72),
226fCalRocArrayOutliers(72),
227fHistoQArray(72),
228fHistoT0Array(72),
229fHistoRMSArray(72),
230fHMeanTimeSector(0x0),
231fVMeanTimeSector(72),
232fPadTimesArrayEvent(72),
233fPadQArrayEvent(72),
234fPadRMSArrayEvent(72),
235fPadPedestalArrayEvent(72),
236fCurrentChannel(-1),
237fCurrentSector(-1),
238fCurrentRow(-1),
239fCurrentPad(-1),
240fMaxPadSignal(-1),
241fMaxTimeBin(-1),
242fPadSignal(1024),
243fPadPedestal(0),
244fPadNoise(0),
245fVTime0Offset(72),
246fVTime0OffsetCounter(72)
467a4b25 247{
880c3382 248 //
249 // AliTPCSignal default constructor
250 //
251 SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
252 fFirstTimeBin=60;
253 fLastTimeBin=900;
ac940b58 254 fParam->Update();
467a4b25 255}
256//_____________________________________________________________________
257AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
c3066940 258AliTPCCalibRawBase(sig),
259fNbinsT0(sig.fNbinsT0),
260fXminT0(sig.fXminT0),
261fXmaxT0(sig.fXmaxT0),
262fNbinsQ(sig.fNbinsQ),
263fXminQ(sig.fXminQ),
264fXmaxQ(sig.fXmaxQ),
265fNbinsRMS(sig.fNbinsRMS),
266fXminRMS(sig.fXminRMS),
267fXmaxRMS(sig.fXmaxRMS),
268fPeakIntMinus(sig.fPeakIntMinus),
269fPeakIntPlus(sig.fPeakIntPlus),
270fIsZeroSuppressed(sig.fIsZeroSuppressed),
271fLastSector(-1),
272fParam(new AliTPCParam),
273fPedestalTPC(0x0),
274fPadNoiseTPC(0x0),
275fOutliers(0x0),
276fPedestalROC(0x0),
277fPadNoiseROC(0x0),
278fCalRocArrayT0(72),
279fCalRocArrayQ(72),
280fCalRocArrayRMS(72),
281fCalRocArrayOutliers(72),
282fHistoQArray(72),
283fHistoT0Array(72),
284fHistoRMSArray(72),
285fHMeanTimeSector(0x0),
286fVMeanTimeSector(72),
287fPadTimesArrayEvent(72),
288fPadQArrayEvent(72),
289fPadRMSArrayEvent(72),
290fPadPedestalArrayEvent(72),
291fCurrentChannel(-1),
292fCurrentSector(-1),
293fCurrentRow(-1),
294fCurrentPad(-1),
295fMaxPadSignal(-1),
296fMaxTimeBin(-1),
297fPadSignal(1024),
298fPadPedestal(0),
299fPadNoise(0),
300fVTime0Offset(72),
301fVTime0OffsetCounter(72)
467a4b25 302{
880c3382 303 //
304 // AliTPCSignal default constructor
305 //
306
307 for (Int_t iSec = 0; iSec < 72; ++iSec){
308 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
309 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
310 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
311 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
312
313 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
314 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
315 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
316
317 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
318 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
319 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
320 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
321
322 if ( hQ != 0x0 ){
323 TH2S *hNew = new TH2S(*hQ);
324 hNew->SetDirectory(0);
325 fHistoQArray.AddAt(hNew,iSec);
467a4b25 326 }
880c3382 327 if ( hT0 != 0x0 ){
328 TH2S *hNew = new TH2S(*hT0);
329 hNew->SetDirectory(0);
330 fHistoQArray.AddAt(hNew,iSec);
331 }
332 if ( hRMS != 0x0 ){
333 TH2S *hNew = new TH2S(*hRMS);
334 hNew->SetDirectory(0);
335 fHistoQArray.AddAt(hNew,iSec);
336 }
337 fVMeanTimeSector[iSec]=sig.fVMeanTimeSector[iSec];
338 }
339
340 if ( sig.fHMeanTimeSector ) fHMeanTimeSector=(TH2F*)sig.fHMeanTimeSector->Clone();
341 fParam->Update();
467a4b25 342}
ac940b58 343AliTPCCalibPulser::AliTPCCalibPulser(const TMap *config) :
c3066940 344AliTPCCalibRawBase(),
345fNbinsT0(200),
346fXminT0(-2),
347fXmaxT0(2),
348fNbinsQ(200),
349fXminQ(10),
350fXmaxQ(40),
351fNbinsRMS(100),
352fXminRMS(0.1),
353fXmaxRMS(5.1),
354fPeakIntMinus(2),
355fPeakIntPlus(2),
356fIsZeroSuppressed(kFALSE),
357fLastSector(-1),
358fParam(new AliTPCParam),
359fPedestalTPC(0x0),
360fPadNoiseTPC(0x0),
361fOutliers(0x0),
362fPedestalROC(0x0),
363fPadNoiseROC(0x0),
364fCalRocArrayT0(72),
365fCalRocArrayQ(72),
366fCalRocArrayRMS(72),
367fCalRocArrayOutliers(72),
368fHistoQArray(72),
369fHistoT0Array(72),
370fHistoRMSArray(72),
371fHMeanTimeSector(0x0),
372fVMeanTimeSector(72),
373fPadTimesArrayEvent(72),
374fPadQArrayEvent(72),
375fPadRMSArrayEvent(72),
376fPadPedestalArrayEvent(72),
377fCurrentChannel(-1),
378fCurrentSector(-1),
379fCurrentRow(-1),
380fCurrentPad(-1),
381fMaxPadSignal(-1),
382fMaxTimeBin(-1),
383fPadSignal(1024),
384fPadPedestal(0),
385fPadNoise(0),
386fVTime0Offset(72),
387fVTime0OffsetCounter(72)
ac940b58 388{
389 //
390 // This constructor uses a TMap for setting some parametes
391 //
880c3382 392 SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
393 fFirstTimeBin=60;
394 fLastTimeBin=900;
ac940b58 395 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
396 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
397 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
398 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
399 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
400 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
401 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
402 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
403 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
404 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
405 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
880c3382 406 if (config->GetValue("PeakIntMinus")) fPeakIntMinus = (Int_t)((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atof();
407 if (config->GetValue("PeakIntPlus")) fPeakIntPlus = (Int_t)((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atof();
ac940b58 408 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
c3066940 409
ac940b58 410 fParam->Update();
411}
467a4b25 412//_____________________________________________________________________
413AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
414{
415 //
416 // assignment operator
417 //
418 if (&source == this) return *this;
419 new (this) AliTPCCalibPulser(source);
c3066940 420
467a4b25 421 return *this;
422}
423//_____________________________________________________________________
424AliTPCCalibPulser::~AliTPCCalibPulser()
425{
880c3382 426 //
427 // destructor
428 //
429
430 Reset();
431 delete fParam;
d8faeea7 432}
433void AliTPCCalibPulser::Reset()
434{
880c3382 435 //
436 // Delete all information: Arrays, Histograms, CalRoc objects
437 //
438 fCalRocArrayT0.Delete();
439 fCalRocArrayQ.Delete();
440 fCalRocArrayRMS.Delete();
441 fCalRocArrayOutliers.Delete();
442
443 fHistoQArray.Delete();
444 fHistoT0Array.Delete();
445 fHistoRMSArray.Delete();
446
447 fPadTimesArrayEvent.Delete();
448 fPadQArrayEvent.Delete();
449 fPadRMSArrayEvent.Delete();
450 fPadPedestalArrayEvent.Delete();
451
ac940b58 452 if (fHMeanTimeSector) delete fHMeanTimeSector;
467a4b25 453}
454//_____________________________________________________________________
3cd27a08 455Int_t AliTPCCalibPulser::Update(const Int_t icsector,
c3066940 456 const Int_t icRow,
457 const Int_t icPad,
458 const Int_t icTimeBin,
459 const Float_t csignal)
467a4b25 460{
461 //
3cd27a08 462 // Signal filling methode on the fly pedestal and time offset correction if necessary.
467a4b25 463 // no extra analysis necessary. Assumes knowledge of the signal shape!
464 // assumes that it is looped over consecutive time bins of one pad
465 //
c3066940 466
467 if (icRow<0) return 0;
468 if (icPad<0) return 0;
469 if (icTimeBin<0) return 0;
470 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
471
472 if ( icRow<0 || icPad<0 ){
473 AliWarning("Wrong Pad or Row number, skipping!");
474 return 0;
475 }
476
477 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
478
467a4b25 479 //init first pad and sector in this event
c3066940 480 if ( fCurrentChannel == -1 ) {
481 fCurrentChannel = iChannel;
482 fCurrentSector = icsector;
483 fCurrentRow = icRow;
484 fCurrentPad = icPad;
485 }
486
467a4b25 487 //process last pad if we change to a new one
c3066940 488 if ( iChannel != fCurrentChannel ){
489 ProcessPad();
490 fLastSector=fCurrentSector;
491 fCurrentChannel = iChannel;
492 fCurrentSector = icsector;
493 fCurrentRow = icRow;
494 fCurrentPad = icPad;
495 }
496
467a4b25 497 //fill signals for current pad
c3066940 498 fPadSignal[icTimeBin]=csignal;
499 if ( csignal > fMaxPadSignal ){
500 fMaxPadSignal = csignal;
501 fMaxTimeBin = icTimeBin;
502 }
503 return 0;
467a4b25 504}
505//_____________________________________________________________________
506void AliTPCCalibPulser::FindPedestal(Float_t part)
507{
508 //
509 // find pedestal and noise for the current pad. Use either database or
510 // truncated mean with part*100%
511 //
c3066940 512 Bool_t noPedestal = kTRUE;;
513 if (fPedestalTPC&&fPadNoiseTPC){
467a4b25 514 //use pedestal database
515 //only load new pedestals if the sector has changed
c3066940 516 if ( fCurrentSector!=fLastSector ){
517 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
518 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
467a4b25 519 }
c3066940 520
521 if ( fPedestalROC&&fPadNoiseROC ){
522 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
523 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
524 noPedestal = kFALSE;
525 }
526
527 }
528
467a4b25 529 //if we are not running with pedestal database, or for the current sector there is no information
530 //available, calculate the pedestal and noise on the fly
c3066940 531 if ( noPedestal ) {
532 const Int_t kPedMax = 100; //maximum pedestal value
533 Float_t max = 0;
c3066940 534 Int_t median = -1;
535 Int_t count0 = 0;
536 Int_t count1 = 0;
537 //
538 Float_t padSignal=0;
467a4b25 539 //
c3066940 540 UShort_t histo[kPedMax];
541 memset(histo,0,kPedMax*sizeof(UShort_t));
542
543 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
544 padSignal = fPadSignal.GetMatrixArray()[i];
545 if (padSignal<=0) continue;
546 if (padSignal>max && i>10) {
547 max = padSignal;
c3066940 548 }
549 if (padSignal>kPedMax-1) continue;
550 histo[Int_t(padSignal+0.5)]++;
551 count0++;
552 }
553 //
554 for (Int_t i=1; i<kPedMax; ++i){
555 if (count1<count0*0.5) median=i;
556 count1+=histo[i];
557 }
558 // truncated mean
559 //
d8faeea7 560 // what if by chance histo[median] == 0 ?!?
c3066940 561 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
562 //
563 for (Int_t idelta=1; idelta<10; ++idelta){
564 if (median-idelta<=0) continue;
565 if (median+idelta>kPedMax) continue;
566 if (count<part*count1){
567 count+=histo[median-idelta];
568 mean +=histo[median-idelta]*(median-idelta);
569 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
570 count+=histo[median+idelta];
571 mean +=histo[median+idelta]*(median+idelta);
572 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
573 }
467a4b25 574 }
c3066940 575 fPadPedestal = 0;
576 fPadNoise = 0;
577 if ( count > 0 ) {
578 mean/=count;
579 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
580 fPadPedestal = mean;
581 fPadNoise = rms;
582 }
583 }
584 fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
467a4b25 585}
586//_____________________________________________________________________
587void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
588{
c3066940 589//
467a4b25 590 // Find position, signal width and height of the CE signal (last signal)
591 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
592 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
593 //
c3066940 594
595 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
596 Int_t cemaxpos = fMaxTimeBin;
597 Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal sum
598 Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal max
599 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
600 const Int_t kCemax = fPeakIntPlus;
601 param[0] = ceQmax;
602 param[1] = ceTime;
603 param[2] = ceRMS;
604 qSum = ceQsum;
605
ac940b58 606 if (cemaxpos>0){
607 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
608 if ( ceQmax<ceMaxThreshold ) return;
880c3382 609 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
ac940b58 610 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
611 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
612 ceTime+=signal*(i+0.5);
613 ceRMS +=signal*(i+0.5)*(i+0.5);
614 ceQsum+=signal;
615 }
467a4b25 616 }
ac940b58 617 }
618 if (ceQsum>ceSumThreshold) {
619 ceTime/=ceQsum;
620 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
880c3382 621 ceTime-=GetL1PhaseTB();
d8faeea7 622 //only fill the Time0Offset if pad was not marked as an outlier!
ac940b58 623 if ( !fOutliers ){
624 //skip edge pads for calculating the mean time
625 if ( !IsEdgePad(fCurrentSector,fCurrentRow,fCurrentPad) ){
626 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
627 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
628 GetHistoTSec()->Fill(ceTime,fCurrentSector);
629 }
630 } else {
631 if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
632 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
633 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
634 }
635 }
636
637 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
638 // the pick-up signal should scale with the pad area. In addition
639 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
640 // ratio 2/3. The pad area we express in mm2 (factor 100). We normalise the signal
641 // to the OROC signal (factor 2/3 for the IROCs).
642 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow)*100;
643 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
c3066940 644
645 ceQsum/=norm;
646 } else {
647 ceQmax=0;
648 ceTime=0;
649 ceRMS =0;
650 ceQsum=0;
651 }
652 param[0] = ceQmax;
653 param[1] = ceTime;
654 param[2] = ceRMS;
655 qSum = ceQsum;
467a4b25 656}
657//_____________________________________________________________________
3cd27a08 658void AliTPCCalibPulser::ProcessPad()
467a4b25 659{
ac940b58 660 //
661 // Process data of current pad
662 //
c3066940 663
ac940b58 664 FindPedestal();
665 TVectorD param(3);
666 Float_t qSum;
667 FindPulserSignal(param, qSum);
668
669 Double_t meanT = param[1];
670 Double_t sigmaT = param[2];
671
672
673 //Fill Event T0 counter
674 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
675
676 //Fill Q histogram
677// GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel ); //use linear scale, needs also a change in the Analyse funciton.
678 GetHistoQ(fCurrentSector,kTRUE)->Fill( qSum, fCurrentChannel );
679
680 //Fill RMS histogram
681 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
682
683
467a4b25 684 //Fill debugging info
880c3382 685 if ( GetStreamLevel()>0 ){
686 TTreeSRedirector *streamer=GetDebugStreamer();
687 if ( GetStreamLevel() == 1 ){
688 if ( streamer ) {
689 Int_t padc = fCurrentPad-(fROC->GetNPads(fCurrentSector,fCurrentRow)/2);
690 (*streamer) << "PadSignals" <<
c3066940 691 "Event=" <<fNevents <<
880c3382 692 "Sector=" <<fCurrentSector<<
693 "Row=" <<fCurrentRow<<
694 "Pad=" <<fCurrentPad<<
695 "PadC=" <<padc<<
696 "Channel="<<fCurrentChannel<<
697 "Sum=" <<qSum<<
698 "params.="<<&param<<
699 "signal.=" <<&fPadSignal<<
700 "\n";
ac940b58 701 }
ac940b58 702 } else { //debug > 1
703 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
704 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
705 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
467a4b25 706 }
ac940b58 707 }
708 ResetPad();
467a4b25 709}
710//_____________________________________________________________________
3cd27a08 711void AliTPCCalibPulser::EndEvent()
467a4b25 712{
713 //
714 // Process data of current event
715 //
ac940b58 716
467a4b25 717 //check if last pad has allready been processed, if not do so
ac940b58 718 if ( fMaxTimeBin>-1 ) ProcessPad();
719
467a4b25 720 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
ac940b58 721 for ( Int_t iSec = 0; iSec<72; ++iSec ){
722 TVectorF *vTimes = GetPadTimesEvent(iSec);
723 if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
724 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
725 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
726 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
727
728 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
d8faeea7 729 //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
ac940b58 730
731
732 //Debug start
880c3382 733 if ( GetStreamLevel()>1 ){
734 TTreeSRedirector *streamer=GetDebugStreamer();
735 if ( streamer ) {
736 Int_t row=0;
737 Int_t pad=0;
738 Int_t padc=0;
739
740 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
741 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
742
743 UInt_t channel=iChannel;
744 Int_t sector=iSec;
745
746 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
747 pad = channel-fROC->GetRowIndexes(sector)[row];
748 padc = pad-(fROC->GetNPads(sector,row)/2);
749
750 (*streamer) << "DataPad" <<
c3066940 751 "Event=" << fNevents <<
880c3382 752 "Sector="<< sector <<
753 "Row=" << row<<
754 "Pad=" << pad <<
755 "PadC=" << padc <<
756 "PadSec="<< channel <<
757 "Time0=" << time0 <<
758 "Time=" << time <<
759 "RMS=" << rms <<
760 "Sum=" << q <<
761 "\n";
ac940b58 762 }
08205ed7 763 }
880c3382 764 //Debug end
4c6d06dc 765 }
08205ed7 766 }
c3066940 767 ++fNevents;
467a4b25 768}
769//_____________________________________________________________________
3cd27a08 770TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
c3066940 771 Int_t nbinsY, Float_t ymin, Float_t ymax,
772 const Char_t *type, Bool_t force)
467a4b25 773{
774 //
775 // return pointer to Q histogram
776 // if force is true create a new histogram if it doesn't exist allready
777 //
c3066940 778 if ( !force || arr->UncheckedAt(sector) )
779 return (TH2S*)arr->UncheckedAt(sector);
780
a3b590cf 781 // if we are forced and histogram doesn't yes exist create it
782 // new histogram with Q calib information. One value for each pad!
783 TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
c3066940 784 nbinsY, ymin, ymax,
785 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
786 hist->SetDirectory(0);
787 arr->AddAt(hist,sector);
788 return hist;
467a4b25 789}
790//_____________________________________________________________________
3cd27a08 791TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
467a4b25 792{
793 //
794 // return pointer to T0 histogram
795 // if force is true create a new histogram if it doesn't exist allready
796 //
c3066940 797 TObjArray *arr = &fHistoT0Array;
798 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
467a4b25 799}
800//_____________________________________________________________________
3cd27a08 801TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
467a4b25 802{
803 //
804 // return pointer to Q histogram
805 // if force is true create a new histogram if it doesn't exist allready
806 //
c3066940 807 TObjArray *arr = &fHistoQArray;
808 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
467a4b25 809}
810//_____________________________________________________________________
3cd27a08 811TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
467a4b25 812{
813 //
814 // return pointer to Q histogram
815 // if force is true create a new histogram if it doesn't exist allready
816 //
c3066940 817 TObjArray *arr = &fHistoRMSArray;
818 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
467a4b25 819}
820//_____________________________________________________________________
ac940b58 821TH2F* AliTPCCalibPulser::GetHistoTSec()
822{
823 //
824 // return the pointer to the abs time distribution per sector
825 // create it if it does not exist
826 //
c3066940 827 if ( !fHMeanTimeSector ) //!!!if you change the binning here, you should also change it in the Analyse Function!!
828 fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
829 20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
830 72,0,72);
831 return fHMeanTimeSector;
ac940b58 832}
833//_____________________________________________________________________
3cd27a08 834TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
467a4b25 835{
836 //
837 // return pointer to Pad Info from 'arr' for the current event and sector
838 // if force is true create it if it doesn't exist allready
839 //
c3066940 840 if ( !force || arr->UncheckedAt(sector) )
841 return (TVectorF*)arr->UncheckedAt(sector);
842
843 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
844 arr->AddAt(vect,sector);
845 return vect;
467a4b25 846}
847//_____________________________________________________________________
3cd27a08 848TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
467a4b25 849{
850 //
851 // return pointer to Pad Times Array for the current event and sector
852 // if force is true create it if it doesn't exist allready
853 //
c3066940 854 TObjArray *arr = &fPadTimesArrayEvent;
855 return GetPadInfoEvent(sector,arr,force);
467a4b25 856}
857//_____________________________________________________________________
3cd27a08 858TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
467a4b25 859{
860 //
861 // return pointer to Pad Q Array for the current event and sector
862 // if force is true create it if it doesn't exist allready
863 // for debugging purposes only
864 //
c3066940 865
866 TObjArray *arr = &fPadQArrayEvent;
867 return GetPadInfoEvent(sector,arr,force);
467a4b25 868}
869//_____________________________________________________________________
3cd27a08 870TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
467a4b25 871{
872 //
873 // return pointer to Pad RMS Array for the current event and sector
874 // if force is true create it if it doesn't exist allready
875 // for debugging purposes only
876 //
c3066940 877 TObjArray *arr = &fPadRMSArrayEvent;
878 return GetPadInfoEvent(sector,arr,force);
467a4b25 879}
880//_____________________________________________________________________
3cd27a08 881TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
467a4b25 882{
883 //
884 // return pointer to Pad RMS Array for the current event and sector
885 // if force is true create it if it doesn't exist allready
886 // for debugging purposes only
887 //
c3066940 888 TObjArray *arr = &fPadPedestalArrayEvent;
889 return GetPadInfoEvent(sector,arr,force);
467a4b25 890}
891//_____________________________________________________________________
3cd27a08 892AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
467a4b25 893{
894 //
895 // return pointer to ROC Calibration
896 // if force is true create a new histogram if it doesn't exist allready
897 //
c3066940 898 if ( !force || arr->UncheckedAt(sector) )
899 return (AliTPCCalROC*)arr->UncheckedAt(sector);
900
467a4b25 901 // if we are forced and histogram doesn't yes exist create it
c3066940 902
467a4b25 903 // new AliTPCCalROC for T0 information. One value for each pad!
c3066940 904 AliTPCCalROC *croc = new AliTPCCalROC(sector);
905 arr->AddAt(croc,sector);
906 return croc;
467a4b25 907}
908//_____________________________________________________________________
3cd27a08 909AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
467a4b25 910{
911 //
912 // return pointer to Carge ROC Calibration
913 // if force is true create a new histogram if it doesn't exist allready
914 //
c3066940 915 TObjArray *arr = &fCalRocArrayT0;
916 return GetCalRoc(sector, arr, force);
467a4b25 917}
918//_____________________________________________________________________
3cd27a08 919AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
467a4b25 920{
921 //
922 // return pointer to T0 ROC Calibration
923 // if force is true create a new histogram if it doesn't exist allready
924 //
c3066940 925 TObjArray *arr = &fCalRocArrayQ;
926 return GetCalRoc(sector, arr, force);
467a4b25 927}
928//_____________________________________________________________________
3cd27a08 929AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
467a4b25 930{
931 //
932 // return pointer to signal width ROC Calibration
933 // if force is true create a new histogram if it doesn't exist allready
934 //
c3066940 935 TObjArray *arr = &fCalRocArrayRMS;
936 return GetCalRoc(sector, arr, force);
467a4b25 937}
938//_____________________________________________________________________
939AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
940{
941 //
942 // return pointer to Outliers
943 // if force is true create a new histogram if it doesn't exist allready
944 //
c3066940 945 TObjArray *arr = &fCalRocArrayOutliers;
946 return GetCalRoc(sector, arr, force);
467a4b25 947}
948//_____________________________________________________________________
3cd27a08 949void AliTPCCalibPulser::ResetEvent()
467a4b25 950{
951 //
952 // Reset global counters -- Should be called before each event is processed
953 //
c3066940 954 fLastSector=-1;
955 fCurrentSector=-1;
956 fCurrentRow=-1;
957 fCurrentPad=-1;
958 fCurrentChannel=-1;
959
960 ResetPad();
961
962 fPadTimesArrayEvent.Delete();
963
964 fPadQArrayEvent.Delete();
965 fPadRMSArrayEvent.Delete();
966 fPadPedestalArrayEvent.Delete();
967
968 for ( Int_t i=0; i<72; ++i ){
969 fVTime0Offset[i]=0;
970 fVTime0OffsetCounter[i]=0;
971 }
467a4b25 972}
973//_____________________________________________________________________
3cd27a08 974void AliTPCCalibPulser::ResetPad()
467a4b25 975{
976 //
977 // Reset pad infos -- Should be called after a pad has been processed
978 //
c3066940 979 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
980 fPadSignal[i] = 0;
981 fMaxTimeBin = -1;
982 fMaxPadSignal = -1;
983 fPadPedestal = -1;
984 fPadNoise = -1;
467a4b25 985}
986//_____________________________________________________________________
ac940b58 987Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
988{
989 //
990 // return true if pad is on the edge of a row
991 //
c3066940 992 Int_t edge1 = 0;
993 Int_t edge2 = fROC->GetNPads(sector,row)-1;
994 if ( pad == edge1 || pad == edge2 ) return kTRUE;
995
996 return kFALSE;
ac940b58 997}
998//_____________________________________________________________________
7442bceb 999void AliTPCCalibPulser::Merge(AliTPCCalibPulser * const sig)
467a4b25 1000{
78f17711 1001 //
1002 // Merge reference histograms of sig to the current AliTPCCalibPulser
1003 //
1004
1005 MergeBase(sig);
1006 //merge histograms
c3066940 1007 for (Int_t iSec=0; iSec<72; ++iSec){
1008 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
1009 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
1010 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1011
1012
1013 if ( hRefQmerge ){
1014 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1015 TH2S *hRefQ = GetHistoQ(iSec);
1016 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1017 else {
1018 TH2S *hist = new TH2S(*hRefQmerge);
1019 hist->SetDirectory(0);
1020 fHistoQArray.AddAt(hist, iSec);
1021 }
1022 hRefQmerge->SetDirectory(dir);
467a4b25 1023 }
c3066940 1024 if ( hRefT0merge ){
1025 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1026 TH2S *hRefT0 = GetHistoT0(iSec);
1027 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1028 else {
1029 TH2S *hist = new TH2S(*hRefT0merge);
1030 hist->SetDirectory(0);
1031 fHistoT0Array.AddAt(hist, iSec);
1032 }
1033 hRefT0merge->SetDirectory(dir);
ac940b58 1034 }
c3066940 1035 if ( hRefRMSmerge ){
1036 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1037 TH2S *hRefRMS = GetHistoRMS(iSec);
1038 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1039 else {
1040 TH2S *hist = new TH2S(*hRefRMSmerge);
1041 hist->SetDirectory(0);
1042 fHistoRMSArray.AddAt(hist, iSec);
1043 }
1044 hRefRMSmerge->SetDirectory(dir);
1045 }
1046
1047 }
1048 if ( sig->fHMeanTimeSector ){
1049 TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1050 if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1051 else {
1052 fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1053 fHMeanTimeSector->SetDirectory(0);
1054 }
1055 sig->fHMeanTimeSector->SetDirectory(dir);
1056 }
467a4b25 1057}
7442bceb 1058
1059
1060//_____________________________________________________________________
1061Long64_t AliTPCCalibPulser::Merge(TCollection * const list)
1062{
1063 //
1064 // Merge all objects of this type in list
1065 //
1066
1067 Long64_t nmerged=1;
1068
1069 TIter next(list);
1070 AliTPCCalibPulser *ce=0;
1071 TObject *o=0;
1072
1073 while ( (o=next()) ){
1074 ce=dynamic_cast<AliTPCCalibPulser*>(o);
1075 if (ce){
1076 Merge(ce);
1077 ++nmerged;
1078 }
1079 }
1080
1081 return nmerged;
1082}
1083
467a4b25 1084//_____________________________________________________________________
3cd27a08 1085void AliTPCCalibPulser::Analyse()
467a4b25 1086{
880c3382 1087 //
1088 // Calculate calibration constants
1089 //
ac940b58 1090
1091 TVectorD paramQ(3);
1092 TVectorD paramT0(3);
1093 TVectorD paramRMS(3);
1094 TMatrixD dummy(3,3);
1095 //calculate mean time for each sector and mean time for each side
1096 TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1097 fVMeanTimeSector.Zero();
1098
1099 for (Int_t iSec=0; iSec<72; ++iSec){
1100 TH2S *hT0 = GetHistoT0(iSec);
1101 if (!hT0 ) continue;
c3066940 1102 //calculate sector mean T
ac940b58 1103 if ( fHMeanTimeSector ){
1104 Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1105 Int_t offset = (nbinsT+2)*(iSec+1);
1106 Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1107 Int_t entries=0;
1108 for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1109 hMeanTsec.Set(nbinsT+2,arrP);
1110 hMeanTsec.SetEntries(entries);
1111 paramT0.Zero();
1112 // truncated mean: remove lower 5% and upper 5%
1113 if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,&paramT0,0.05,.95);
1114 fVMeanTimeSector[iSec]=paramT0[1];
c3066940 1115 }
ac940b58 1116
1117 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1118 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1119 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1120 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1121
1122 TH2S *hQ = GetHistoQ(iSec);
1123 TH2S *hRMS = GetHistoRMS(iSec);
1124
1125 Short_t *arrayhQ = hQ->GetArray();
1126 Short_t *arrayhT0 = hT0->GetArray();
1127 Short_t *arrayhRMS = hRMS->GetArray();
1128
1129 UInt_t nChannels = fROC->GetNChannels(iSec);
1130 Float_t meanTsec = fVMeanTimeSector[iSec];
c3066940 1131
ac940b58 1132 //debug
1133 Int_t row=0;
1134 Int_t pad=0;
1135 Int_t padc=0;
1136 //! debug
1137
1138 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1139
1140 Float_t cogTime0 = -1000;
1141 Float_t cogQ = -1000;
1142 Float_t cogRMS = -1000;
1143 Float_t cogOut = 0;
1144
1145 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1146 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1147 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
467a4b25 1148/*
ac940b58 1149 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1150 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1151 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1152 cogQ = paramQ[1];
1153 cogTime0 = paramT0[1];
1154 cogRMS = paramRMS[1];
467a4b25 1155*/
ac940b58 1156 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1157 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1158 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1159
1160 /*
1161 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1162 cogOut = 1;
1163 cogTime0 = 0;
1164 cogQ = 0;
1165 cogRMS = 0;
1166 }
467a4b25 1167*/
ac940b58 1168// rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1169 rocQ->SetValue(iChannel, cogQ);
1170 rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1171 rocRMS->SetValue(iChannel, cogRMS);
1172 rocOut->SetValue(iChannel, cogOut);
12f70d37 1173
1174 // in case a channel has no data set the value to 0
1175 if (TMath::Abs(cogTime0-fXminT0)<1e-10){
1176 rocQ->SetValue(iChannel, 0);
1177 rocT0->SetValue(iChannel, 0); //offset by mean time of the sector
1178 rocRMS->SetValue(iChannel, 0);
1179 }
1180
ac940b58 1181 //debug
880c3382 1182 if ( GetStreamLevel() > 2 ){
1183 TTreeSRedirector *streamer=GetDebugStreamer();
1184 if ( streamer ) {
1185 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1186 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1187 padc = pad-(fROC->GetNPads(iSec,row)/2);
1188
1189 (*streamer) << "DataEnd" <<
1190 "Sector=" << iSec <<
1191 "Pad=" << pad <<
1192 "PadC=" << padc <<
1193 "Row=" << row <<
1194 "PadSec=" << iChannel <<
1195 "Q=" << cogQ <<
1196 "T0=" << cogTime0 <<
1197 "RMS=" << cogRMS <<
1198 "\n";
ac940b58 1199 }
ac940b58 1200 }
1201 //! debug
467a4b25 1202 }
ac940b58 1203
1204
1205 }
467a4b25 1206}
1207//_____________________________________________________________________
1208//_________________________ Test Functions ___________________________
1209//_____________________________________________________________________
1210TObjArray* AliTPCCalibPulser::TestBinning()
1211{
ac940b58 1212 //
1213 // Function to test the binning of the reference histograms
1214 // type: T0, Q or RMS
1215 // mode: 0 - number of filled bins per channel
1216 // 1 - number of empty bins between filled bins in one ROC
1217 // returns TObjArray with the test histograms type*2+mode:
1218 // position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1219
1220
1221 TObjArray *histArray = new TObjArray(6);
1222 const Char_t *type[] = {"T0","Q","RMS"};
1223 Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1224
1225 for (Int_t itype = 0; itype<3; ++itype){
1226 for (Int_t imode=0; imode<2; ++imode){
1227 Int_t icount = itype*2+imode;
1228 histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1229 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1230 72,0,72),
1231 icount);
467a4b25 1232 }
ac940b58 1233 }
1234
1235
1236 TH2S *hRef=0x0;
1237 Short_t *array=0x0;
1238 for (Int_t itype = 0; itype<3; ++itype){
1239 for (Int_t iSec=0; iSec<72; ++iSec){
1240 if ( itype == 0 ) hRef = GetHistoT0(iSec);
1241 if ( itype == 1 ) hRef = GetHistoQ(iSec);
1242 if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1243 if ( hRef == 0x0 ) continue;
1244 array = (hRef->GetArray());
1245 UInt_t nChannels = fROC->GetNChannels(iSec);
1246
1247 Int_t nempty=0;
1248 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1249 Int_t nfilled=0;
1250 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1251 Int_t c1 = 0;
1252 Int_t c2 = 0;
1253 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1254 if ( array[offset+iBin]>0 ) {
1255 nfilled++;
1256 if ( c1 && c2 ) nempty++;
1257 else c1 = 1;
1258 }
1259 else if ( c1 ) c2 = 1;
1260
1261 }
1262 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1263 }
1264 ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
467a4b25 1265 }
ac940b58 1266 }
1267 return histArray;
467a4b25 1268}