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