15 minutes interval for calculation of drift correction
[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() :
c3066940 205AliTPCCalibRawBase(),
206fNbinsT0(200),
207fXminT0(-2),
208fXmaxT0(2),
209fNbinsQ(200),
210fXminQ(10),
211fXmaxQ(40),
212fNbinsRMS(100),
213fXminRMS(0.1),
214fXmaxRMS(5.1),
215fPeakIntMinus(2),
216fPeakIntPlus(2),
217fIsZeroSuppressed(kFALSE),
218fLastSector(-1),
219fParam(new AliTPCParam),
220fPedestalTPC(0x0),
221fPadNoiseTPC(0x0),
222fOutliers(0x0),
223fPedestalROC(0x0),
224fPadNoiseROC(0x0),
225fCalRocArrayT0(72),
226fCalRocArrayQ(72),
227fCalRocArrayRMS(72),
228fCalRocArrayOutliers(72),
229fHistoQArray(72),
230fHistoT0Array(72),
231fHistoRMSArray(72),
232fHMeanTimeSector(0x0),
233fVMeanTimeSector(72),
234fPadTimesArrayEvent(72),
235fPadQArrayEvent(72),
236fPadRMSArrayEvent(72),
237fPadPedestalArrayEvent(72),
238fCurrentChannel(-1),
239fCurrentSector(-1),
240fCurrentRow(-1),
241fCurrentPad(-1),
242fMaxPadSignal(-1),
243fMaxTimeBin(-1),
244fPadSignal(1024),
245fPadPedestal(0),
246fPadNoise(0),
247fVTime0Offset(72),
248fVTime0OffsetCounter(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) :
c3066940 260AliTPCCalibRawBase(sig),
261fNbinsT0(sig.fNbinsT0),
262fXminT0(sig.fXminT0),
263fXmaxT0(sig.fXmaxT0),
264fNbinsQ(sig.fNbinsQ),
265fXminQ(sig.fXminQ),
266fXmaxQ(sig.fXmaxQ),
267fNbinsRMS(sig.fNbinsRMS),
268fXminRMS(sig.fXminRMS),
269fXmaxRMS(sig.fXmaxRMS),
270fPeakIntMinus(sig.fPeakIntMinus),
271fPeakIntPlus(sig.fPeakIntPlus),
272fIsZeroSuppressed(sig.fIsZeroSuppressed),
273fLastSector(-1),
274fParam(new AliTPCParam),
275fPedestalTPC(0x0),
276fPadNoiseTPC(0x0),
277fOutliers(0x0),
278fPedestalROC(0x0),
279fPadNoiseROC(0x0),
280fCalRocArrayT0(72),
281fCalRocArrayQ(72),
282fCalRocArrayRMS(72),
283fCalRocArrayOutliers(72),
284fHistoQArray(72),
285fHistoT0Array(72),
286fHistoRMSArray(72),
287fHMeanTimeSector(0x0),
288fVMeanTimeSector(72),
289fPadTimesArrayEvent(72),
290fPadQArrayEvent(72),
291fPadRMSArrayEvent(72),
292fPadPedestalArrayEvent(72),
293fCurrentChannel(-1),
294fCurrentSector(-1),
295fCurrentRow(-1),
296fCurrentPad(-1),
297fMaxPadSignal(-1),
298fMaxTimeBin(-1),
299fPadSignal(1024),
300fPadPedestal(0),
301fPadNoise(0),
302fVTime0Offset(72),
303fVTime0OffsetCounter(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) :
c3066940 346AliTPCCalibRawBase(),
347fNbinsT0(200),
348fXminT0(-2),
349fXmaxT0(2),
350fNbinsQ(200),
351fXminQ(10),
352fXmaxQ(40),
353fNbinsRMS(100),
354fXminRMS(0.1),
355fXmaxRMS(5.1),
356fPeakIntMinus(2),
357fPeakIntPlus(2),
358fIsZeroSuppressed(kFALSE),
359fLastSector(-1),
360fParam(new AliTPCParam),
361fPedestalTPC(0x0),
362fPadNoiseTPC(0x0),
363fOutliers(0x0),
364fPedestalROC(0x0),
365fPadNoiseROC(0x0),
366fCalRocArrayT0(72),
367fCalRocArrayQ(72),
368fCalRocArrayRMS(72),
369fCalRocArrayOutliers(72),
370fHistoQArray(72),
371fHistoT0Array(72),
372fHistoRMSArray(72),
373fHMeanTimeSector(0x0),
374fVMeanTimeSector(72),
375fPadTimesArrayEvent(72),
376fPadQArrayEvent(72),
377fPadRMSArrayEvent(72),
378fPadPedestalArrayEvent(72),
379fCurrentChannel(-1),
380fCurrentSector(-1),
381fCurrentRow(-1),
382fCurrentPad(-1),
383fMaxPadSignal(-1),
384fMaxTimeBin(-1),
385fPadSignal(1024),
386fPadPedestal(0),
387fPadNoise(0),
388fVTime0Offset(72),
389fVTime0OffsetCounter(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();
c3066940 411
ac940b58 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);
c3066940 422
467a4b25 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,
c3066940 458 const Int_t icRow,
459 const Int_t icPad,
460 const Int_t icTimeBin,
461 const Float_t csignal)
467a4b25 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 //
c3066940 468
469 if (icRow<0) return 0;
470 if (icPad<0) return 0;
471 if (icTimeBin<0) return 0;
472 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
473
474 if ( icRow<0 || icPad<0 ){
475 AliWarning("Wrong Pad or Row number, skipping!");
476 return 0;
477 }
478
479 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
480
467a4b25 481 //init first pad and sector in this event
c3066940 482 if ( fCurrentChannel == -1 ) {
483 fCurrentChannel = iChannel;
484 fCurrentSector = icsector;
485 fCurrentRow = icRow;
486 fCurrentPad = icPad;
487 }
488
467a4b25 489 //process last pad if we change to a new one
c3066940 490 if ( iChannel != fCurrentChannel ){
491 ProcessPad();
492 fLastSector=fCurrentSector;
493 fCurrentChannel = iChannel;
494 fCurrentSector = icsector;
495 fCurrentRow = icRow;
496 fCurrentPad = icPad;
497 }
498
467a4b25 499 //fill signals for current pad
c3066940 500 fPadSignal[icTimeBin]=csignal;
501 if ( csignal > fMaxPadSignal ){
502 fMaxPadSignal = csignal;
503 fMaxTimeBin = icTimeBin;
504 }
505 return 0;
467a4b25 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 //
c3066940 514 Bool_t noPedestal = kTRUE;;
515 if (fPedestalTPC&&fPadNoiseTPC){
467a4b25 516 //use pedestal database
517 //only load new pedestals if the sector has changed
c3066940 518 if ( fCurrentSector!=fLastSector ){
519 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
520 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
467a4b25 521 }
c3066940 522
523 if ( fPedestalROC&&fPadNoiseROC ){
524 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
525 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
526 noPedestal = kFALSE;
527 }
528
529 }
530
467a4b25 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
c3066940 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;
467a4b25 542 //
c3066940 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 ?!?
c3066940 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 }
467a4b25 578 }
c3066940 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 }
588 fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
467a4b25 589}
590//_____________________________________________________________________
591void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
592{
c3066940 593//
467a4b25 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 //
c3066940 598
599 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
600 Int_t cemaxpos = fMaxTimeBin;
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
603 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
604 const Int_t kCemax = fPeakIntPlus;
605 param[0] = ceQmax;
606 param[1] = ceTime;
607 param[2] = ceRMS;
608 qSum = ceQsum;
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.;
c3066940 648
649 ceQsum/=norm;
650 } else {
651 ceQmax=0;
652 ceTime=0;
653 ceRMS =0;
654 ceQsum=0;
655 }
656 param[0] = ceQmax;
657 param[1] = ceTime;
658 param[2] = ceRMS;
659 qSum = ceQsum;
467a4b25 660}
661//_____________________________________________________________________
3cd27a08 662void AliTPCCalibPulser::ProcessPad()
467a4b25 663{
ac940b58 664 //
665 // Process data of current pad
666 //
c3066940 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" <<
c3066940 695 "Event=" <<fNevents <<
880c3382 696 "Sector=" <<fCurrentSector<<
697 "Row=" <<fCurrentRow<<
698 "Pad=" <<fCurrentPad<<
699 "PadC=" <<padc<<
700 "Channel="<<fCurrentChannel<<
701 "Sum=" <<qSum<<
702 "params.="<<&param<<
703 "signal.=" <<&fPadSignal<<
704 "\n";
ac940b58 705 }
ac940b58 706 } else { //debug > 1
707 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
708 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
709 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
467a4b25 710 }
ac940b58 711 }
712 ResetPad();
467a4b25 713}
714//_____________________________________________________________________
3cd27a08 715void AliTPCCalibPulser::EndEvent()
467a4b25 716{
717 //
718 // Process data of current event
719 //
ac940b58 720
467a4b25 721 //check if last pad has allready been processed, if not do so
ac940b58 722 if ( fMaxTimeBin>-1 ) ProcessPad();
723
467a4b25 724 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
ac940b58 725 for ( Int_t iSec = 0; iSec<72; ++iSec ){
726 TVectorF *vTimes = GetPadTimesEvent(iSec);
727 if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
728 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
729 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
730 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
731
732 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
d8faeea7 733 //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
ac940b58 734
735
736 //Debug start
880c3382 737 if ( GetStreamLevel()>1 ){
738 TTreeSRedirector *streamer=GetDebugStreamer();
739 if ( streamer ) {
740 Int_t row=0;
741 Int_t pad=0;
742 Int_t padc=0;
743
744 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
745 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
746
747 UInt_t channel=iChannel;
748 Int_t sector=iSec;
749
750 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
751 pad = channel-fROC->GetRowIndexes(sector)[row];
752 padc = pad-(fROC->GetNPads(sector,row)/2);
753
754 (*streamer) << "DataPad" <<
c3066940 755 "Event=" << fNevents <<
880c3382 756 "Sector="<< sector <<
757 "Row=" << row<<
758 "Pad=" << pad <<
759 "PadC=" << padc <<
760 "PadSec="<< channel <<
761 "Time0=" << time0 <<
762 "Time=" << time <<
763 "RMS=" << rms <<
764 "Sum=" << q <<
765 "\n";
ac940b58 766 }
08205ed7 767 }
880c3382 768 //Debug end
4c6d06dc 769 }
08205ed7 770 }
c3066940 771 ++fNevents;
467a4b25 772}
773//_____________________________________________________________________
3cd27a08 774TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
c3066940 775 Int_t nbinsY, Float_t ymin, Float_t ymax,
776 const Char_t *type, Bool_t force)
467a4b25 777{
778 //
779 // return pointer to Q histogram
780 // if force is true create a new histogram if it doesn't exist allready
781 //
c3066940 782 if ( !force || arr->UncheckedAt(sector) )
783 return (TH2S*)arr->UncheckedAt(sector);
784
467a4b25 785 // if we are forced and histogram doesn't yes exist create it
c3066940 786 Char_t name[255], title[255];
787
788 sprintf(name,"hCalib%s%.2d",type,sector);
789 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
790
467a4b25 791 // new histogram with Q calib information. One value for each pad!
c3066940 792 TH2S* hist = new TH2S(name,title,
793 nbinsY, ymin, ymax,
794 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
795 hist->SetDirectory(0);
796 arr->AddAt(hist,sector);
797 return hist;
467a4b25 798}
799//_____________________________________________________________________
3cd27a08 800TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
467a4b25 801{
802 //
803 // return pointer to T0 histogram
804 // if force is true create a new histogram if it doesn't exist allready
805 //
c3066940 806 TObjArray *arr = &fHistoT0Array;
807 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
467a4b25 808}
809//_____________________________________________________________________
3cd27a08 810TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
467a4b25 811{
812 //
813 // return pointer to Q histogram
814 // if force is true create a new histogram if it doesn't exist allready
815 //
c3066940 816 TObjArray *arr = &fHistoQArray;
817 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
467a4b25 818}
819//_____________________________________________________________________
3cd27a08 820TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
467a4b25 821{
822 //
823 // return pointer to Q histogram
824 // if force is true create a new histogram if it doesn't exist allready
825 //
c3066940 826 TObjArray *arr = &fHistoRMSArray;
827 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
467a4b25 828}
829//_____________________________________________________________________
ac940b58 830TH2F* AliTPCCalibPulser::GetHistoTSec()
831{
832 //
833 // return the pointer to the abs time distribution per sector
834 // create it if it does not exist
835 //
c3066940 836 if ( !fHMeanTimeSector ) //!!!if you change the binning here, you should also change it in the Analyse Function!!
837 fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
838 20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
839 72,0,72);
840 return fHMeanTimeSector;
ac940b58 841}
842//_____________________________________________________________________
3cd27a08 843TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
467a4b25 844{
845 //
846 // return pointer to Pad Info from 'arr' for the current event and sector
847 // if force is true create it if it doesn't exist allready
848 //
c3066940 849 if ( !force || arr->UncheckedAt(sector) )
850 return (TVectorF*)arr->UncheckedAt(sector);
851
852 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
853 arr->AddAt(vect,sector);
854 return vect;
467a4b25 855}
856//_____________________________________________________________________
3cd27a08 857TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
467a4b25 858{
859 //
860 // return pointer to Pad Times Array for the current event and sector
861 // if force is true create it if it doesn't exist allready
862 //
c3066940 863 TObjArray *arr = &fPadTimesArrayEvent;
864 return GetPadInfoEvent(sector,arr,force);
467a4b25 865}
866//_____________________________________________________________________
3cd27a08 867TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
467a4b25 868{
869 //
870 // return pointer to Pad Q Array for the current event and sector
871 // if force is true create it if it doesn't exist allready
872 // for debugging purposes only
873 //
c3066940 874
875 TObjArray *arr = &fPadQArrayEvent;
876 return GetPadInfoEvent(sector,arr,force);
467a4b25 877}
878//_____________________________________________________________________
3cd27a08 879TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
467a4b25 880{
881 //
882 // return pointer to Pad RMS Array for the current event and sector
883 // if force is true create it if it doesn't exist allready
884 // for debugging purposes only
885 //
c3066940 886 TObjArray *arr = &fPadRMSArrayEvent;
887 return GetPadInfoEvent(sector,arr,force);
467a4b25 888}
889//_____________________________________________________________________
3cd27a08 890TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
467a4b25 891{
892 //
893 // return pointer to Pad RMS Array for the current event and sector
894 // if force is true create it if it doesn't exist allready
895 // for debugging purposes only
896 //
c3066940 897 TObjArray *arr = &fPadPedestalArrayEvent;
898 return GetPadInfoEvent(sector,arr,force);
467a4b25 899}
900//_____________________________________________________________________
3cd27a08 901AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
467a4b25 902{
903 //
904 // return pointer to ROC Calibration
905 // if force is true create a new histogram if it doesn't exist allready
906 //
c3066940 907 if ( !force || arr->UncheckedAt(sector) )
908 return (AliTPCCalROC*)arr->UncheckedAt(sector);
909
467a4b25 910 // if we are forced and histogram doesn't yes exist create it
c3066940 911
467a4b25 912 // new AliTPCCalROC for T0 information. One value for each pad!
c3066940 913 AliTPCCalROC *croc = new AliTPCCalROC(sector);
914 arr->AddAt(croc,sector);
915 return croc;
467a4b25 916}
917//_____________________________________________________________________
3cd27a08 918AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
467a4b25 919{
920 //
921 // return pointer to Carge ROC Calibration
922 // if force is true create a new histogram if it doesn't exist allready
923 //
c3066940 924 TObjArray *arr = &fCalRocArrayT0;
925 return GetCalRoc(sector, arr, force);
467a4b25 926}
927//_____________________________________________________________________
3cd27a08 928AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
467a4b25 929{
930 //
931 // return pointer to T0 ROC Calibration
932 // if force is true create a new histogram if it doesn't exist allready
933 //
c3066940 934 TObjArray *arr = &fCalRocArrayQ;
935 return GetCalRoc(sector, arr, force);
467a4b25 936}
937//_____________________________________________________________________
3cd27a08 938AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
467a4b25 939{
940 //
941 // return pointer to signal width ROC Calibration
942 // if force is true create a new histogram if it doesn't exist allready
943 //
c3066940 944 TObjArray *arr = &fCalRocArrayRMS;
945 return GetCalRoc(sector, arr, force);
467a4b25 946}
947//_____________________________________________________________________
948AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
949{
950 //
951 // return pointer to Outliers
952 // if force is true create a new histogram if it doesn't exist allready
953 //
c3066940 954 TObjArray *arr = &fCalRocArrayOutliers;
955 return GetCalRoc(sector, arr, force);
467a4b25 956}
957//_____________________________________________________________________
3cd27a08 958void AliTPCCalibPulser::ResetEvent()
467a4b25 959{
960 //
961 // Reset global counters -- Should be called before each event is processed
962 //
c3066940 963 fLastSector=-1;
964 fCurrentSector=-1;
965 fCurrentRow=-1;
966 fCurrentPad=-1;
967 fCurrentChannel=-1;
968
969 ResetPad();
970
971 fPadTimesArrayEvent.Delete();
972
973 fPadQArrayEvent.Delete();
974 fPadRMSArrayEvent.Delete();
975 fPadPedestalArrayEvent.Delete();
976
977 for ( Int_t i=0; i<72; ++i ){
978 fVTime0Offset[i]=0;
979 fVTime0OffsetCounter[i]=0;
980 }
467a4b25 981}
982//_____________________________________________________________________
3cd27a08 983void AliTPCCalibPulser::ResetPad()
467a4b25 984{
985 //
986 // Reset pad infos -- Should be called after a pad has been processed
987 //
c3066940 988 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
989 fPadSignal[i] = 0;
990 fMaxTimeBin = -1;
991 fMaxPadSignal = -1;
992 fPadPedestal = -1;
993 fPadNoise = -1;
467a4b25 994}
995//_____________________________________________________________________
ac940b58 996Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
997{
998 //
999 // return true if pad is on the edge of a row
1000 //
c3066940 1001 Int_t edge1 = 0;
1002 Int_t edge2 = fROC->GetNPads(sector,row)-1;
1003 if ( pad == edge1 || pad == edge2 ) return kTRUE;
1004
1005 return kFALSE;
ac940b58 1006}
1007//_____________________________________________________________________
7442bceb 1008void AliTPCCalibPulser::Merge(AliTPCCalibPulser * const sig)
467a4b25 1009{
1010 //
1011 // Merge reference histograms of sig to the current AliTPCCalibPulser
1012 //
c3066940 1013
467a4b25 1014 //merge histograms
c3066940 1015 for (Int_t iSec=0; iSec<72; ++iSec){
1016 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
1017 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
1018 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1019
1020
1021 if ( hRefQmerge ){
1022 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1023 TH2S *hRefQ = GetHistoQ(iSec);
1024 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1025 else {
1026 TH2S *hist = new TH2S(*hRefQmerge);
1027 hist->SetDirectory(0);
1028 fHistoQArray.AddAt(hist, iSec);
1029 }
1030 hRefQmerge->SetDirectory(dir);
467a4b25 1031 }
c3066940 1032 if ( hRefT0merge ){
1033 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1034 TH2S *hRefT0 = GetHistoT0(iSec);
1035 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1036 else {
1037 TH2S *hist = new TH2S(*hRefT0merge);
1038 hist->SetDirectory(0);
1039 fHistoT0Array.AddAt(hist, iSec);
1040 }
1041 hRefT0merge->SetDirectory(dir);
ac940b58 1042 }
c3066940 1043 if ( hRefRMSmerge ){
1044 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1045 TH2S *hRefRMS = GetHistoRMS(iSec);
1046 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1047 else {
1048 TH2S *hist = new TH2S(*hRefRMSmerge);
1049 hist->SetDirectory(0);
1050 fHistoRMSArray.AddAt(hist, iSec);
1051 }
1052 hRefRMSmerge->SetDirectory(dir);
1053 }
1054
1055 }
1056 if ( sig->fHMeanTimeSector ){
1057 TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1058 if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1059 else {
1060 fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1061 fHMeanTimeSector->SetDirectory(0);
1062 }
1063 sig->fHMeanTimeSector->SetDirectory(dir);
1064 }
467a4b25 1065}
7442bceb 1066
1067
1068//_____________________________________________________________________
1069Long64_t AliTPCCalibPulser::Merge(TCollection * const list)
1070{
1071 //
1072 // Merge all objects of this type in list
1073 //
1074
1075 Long64_t nmerged=1;
1076
1077 TIter next(list);
1078 AliTPCCalibPulser *ce=0;
1079 TObject *o=0;
1080
1081 while ( (o=next()) ){
1082 ce=dynamic_cast<AliTPCCalibPulser*>(o);
1083 if (ce){
1084 Merge(ce);
1085 ++nmerged;
1086 }
1087 }
1088
1089 return nmerged;
1090}
1091
467a4b25 1092//_____________________________________________________________________
3cd27a08 1093void AliTPCCalibPulser::Analyse()
467a4b25 1094{
880c3382 1095 //
1096 // Calculate calibration constants
1097 //
ac940b58 1098
1099 TVectorD paramQ(3);
1100 TVectorD paramT0(3);
1101 TVectorD paramRMS(3);
1102 TMatrixD dummy(3,3);
1103 //calculate mean time for each sector and mean time for each side
1104 TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1105 fVMeanTimeSector.Zero();
1106
1107 for (Int_t iSec=0; iSec<72; ++iSec){
1108 TH2S *hT0 = GetHistoT0(iSec);
1109 if (!hT0 ) continue;
c3066940 1110 //calculate sector mean T
ac940b58 1111 if ( fHMeanTimeSector ){
1112 Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1113 Int_t offset = (nbinsT+2)*(iSec+1);
1114 Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1115 Int_t entries=0;
1116 for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1117 hMeanTsec.Set(nbinsT+2,arrP);
1118 hMeanTsec.SetEntries(entries);
1119 paramT0.Zero();
1120 // truncated mean: remove lower 5% and upper 5%
1121 if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,&paramT0,0.05,.95);
1122 fVMeanTimeSector[iSec]=paramT0[1];
c3066940 1123 }
ac940b58 1124
1125 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1126 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1127 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1128 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1129
1130 TH2S *hQ = GetHistoQ(iSec);
1131 TH2S *hRMS = GetHistoRMS(iSec);
1132
1133 Short_t *arrayhQ = hQ->GetArray();
1134 Short_t *arrayhT0 = hT0->GetArray();
1135 Short_t *arrayhRMS = hRMS->GetArray();
1136
1137 UInt_t nChannels = fROC->GetNChannels(iSec);
1138 Float_t meanTsec = fVMeanTimeSector[iSec];
c3066940 1139
ac940b58 1140 //debug
1141 Int_t row=0;
1142 Int_t pad=0;
1143 Int_t padc=0;
1144 //! debug
1145
1146 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1147
1148 Float_t cogTime0 = -1000;
1149 Float_t cogQ = -1000;
1150 Float_t cogRMS = -1000;
1151 Float_t cogOut = 0;
1152
1153 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1154 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1155 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
467a4b25 1156/*
ac940b58 1157 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1158 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1159 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1160 cogQ = paramQ[1];
1161 cogTime0 = paramT0[1];
1162 cogRMS = paramRMS[1];
467a4b25 1163*/
ac940b58 1164 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1165 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1166 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1167
1168 /*
1169 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1170 cogOut = 1;
1171 cogTime0 = 0;
1172 cogQ = 0;
1173 cogRMS = 0;
1174 }
467a4b25 1175*/
ac940b58 1176// rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1177 rocQ->SetValue(iChannel, cogQ);
1178 rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1179 rocRMS->SetValue(iChannel, cogRMS);
1180 rocOut->SetValue(iChannel, cogOut);
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}