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