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