]> 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
Using the Altro mapping from calib DB
[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>
169#include <TH2S.h>
170#include <TString.h>
171#include <TVectorF.h>
172#include <TMath.h>
173
174#include <TDirectory.h>
175#include <TSystem.h>
176#include <TROOT.h>
177#include <TFile.h>
178
179//AliRoot includes
180#include "AliRawReader.h"
181#include "AliRawReaderRoot.h"
182#include "AliRawReaderDate.h"
183#include "AliTPCRawStream.h"
184#include "AliTPCCalROC.h"
185#include "AliTPCCalPad.h"
186#include "AliTPCROC.h"
187#include "AliTPCParam.h"
188#include "AliTPCCalibPulser.h"
189#include "AliTPCcalibDB.h"
190#include "AliMathBase.h"
191#include "TTreeStream.h"
467a4b25 192
3cd27a08 193//date
194#include "event.h"
195
196
197
198
199ClassImp(AliTPCCalibPulser)
200
201AliTPCCalibPulser::AliTPCCalibPulser() :
467a4b25 202 TObject(),
203 fFirstTimeBin(60),
204 fLastTimeBin(120),
205 fNbinsT0(200),
206 fXminT0(-2),
207 fXmaxT0(2),
208 fNbinsQ(200),
209 fXminQ(1),
210 fXmaxQ(40),
211 fNbinsRMS(100),
212 fXminRMS(0.1),
213 fXmaxRMS(5.1),
214 fLastSector(-1),
215 fOldRCUformat(kTRUE),
216 fROC(AliTPCROC::Instance()),
217 fParam(new AliTPCParam),
218 fPedestalTPC(0x0),
219 fPadNoiseTPC(0x0),
220 fPedestalROC(0x0),
221 fPadNoiseROC(0x0),
222 fCalRocArrayT0(72),
223 fCalRocArrayQ(72),
224 fCalRocArrayRMS(72),
225 fCalRocArrayOutliers(72),
226 fHistoQArray(72),
227 fHistoT0Array(72),
228 fHistoRMSArray(72),
229 fPadTimesArrayEvent(72),
230 fPadQArrayEvent(72),
231 fPadRMSArrayEvent(72),
232 fPadPedestalArrayEvent(72),
233 fCurrentChannel(-1),
234 fCurrentSector(-1),
235 fCurrentRow(-1),
236 fMaxPadSignal(-1),
237 fMaxTimeBin(-1),
238 fPadSignal(1024),
239 fPadPedestal(0),
240 fPadNoise(0),
241 fVTime0Offset(72),
242 fVTime0OffsetCounter(72),
3cd27a08 243// fEvent(-1),
467a4b25 244 fDebugStreamer(0x0),
245 fDebugLevel(0)
246{
247 //
248 // AliTPCSignal default constructor
249 //
250
251}
252//_____________________________________________________________________
253AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
254 TObject(sig),
255 fFirstTimeBin(sig.fFirstTimeBin),
256 fLastTimeBin(sig.fLastTimeBin),
257 fNbinsT0(sig.fNbinsT0),
258 fXminT0(sig.fXminT0),
259 fXmaxT0(sig.fXmaxT0),
260 fNbinsQ(sig.fNbinsQ),
261 fXminQ(sig.fXminQ),
262 fXmaxQ(sig.fXmaxQ),
263 fNbinsRMS(sig.fNbinsRMS),
264 fXminRMS(sig.fXminRMS),
265 fXmaxRMS(sig.fXmaxRMS),
266 fLastSector(-1),
267 fOldRCUformat(kTRUE),
268 fROC(AliTPCROC::Instance()),
269 fParam(new AliTPCParam),
270 fPedestalTPC(0x0),
271 fPadNoiseTPC(0x0),
272 fPedestalROC(0x0),
273 fPadNoiseROC(0x0),
274 fCalRocArrayT0(72),
275 fCalRocArrayQ(72),
276 fCalRocArrayRMS(72),
277 fCalRocArrayOutliers(72),
278 fHistoQArray(72),
279 fHistoT0Array(72),
280 fHistoRMSArray(72),
281 fPadTimesArrayEvent(72),
282 fPadQArrayEvent(72),
283 fPadRMSArrayEvent(72),
284 fPadPedestalArrayEvent(72),
285 fCurrentChannel(-1),
286 fCurrentSector(-1),
287 fCurrentRow(-1),
288 fMaxPadSignal(-1),
289 fMaxTimeBin(-1),
290 fPadSignal(1024),
291 fPadPedestal(0),
292 fPadNoise(0),
293 fVTime0Offset(72),
294 fVTime0OffsetCounter(72),
3cd27a08 295// fEvent(-1),
467a4b25 296 fDebugStreamer(0x0),
297 fDebugLevel(sig.fDebugLevel)
298{
299 //
300 // AliTPCSignal default constructor
301 //
302
303 for (Int_t iSec = 0; iSec < 72; ++iSec){
304 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
305 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
306 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
307 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
308
309 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
310 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
311 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
312
313 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
314 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
315 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
316 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
317
318 if ( hQ != 0x0 ){
319 TH2S *hNew = new TH2S(*hQ);
320 hNew->SetDirectory(0);
321 fHistoQArray.AddAt(hNew,iSec);
322 }
323 if ( hT0 != 0x0 ){
324 TH2S *hNew = new TH2S(*hT0);
325 hNew->SetDirectory(0);
326 fHistoQArray.AddAt(hNew,iSec);
327 }
328 if ( hRMS != 0x0 ){
329 TH2S *hNew = new TH2S(*hRMS);
330 hNew->SetDirectory(0);
331 fHistoQArray.AddAt(hNew,iSec);
332 }
333 }
334
335}
336//_____________________________________________________________________
337AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
338{
339 //
340 // assignment operator
341 //
342 if (&source == this) return *this;
343 new (this) AliTPCCalibPulser(source);
344
345 return *this;
346}
347//_____________________________________________________________________
348AliTPCCalibPulser::~AliTPCCalibPulser()
349{
350 //
351 // destructor
352 //
353
354 fCalRocArrayT0.Delete();
355 fCalRocArrayQ.Delete();
356 fCalRocArrayRMS.Delete();
357 fCalRocArrayOutliers.Delete();
358
359 fHistoQArray.Delete();
360 fHistoT0Array.Delete();
361 fHistoRMSArray.Delete();
362
363 fPadTimesArrayEvent.Delete();
364 fPadQArrayEvent.Delete();
365 fPadRMSArrayEvent.Delete();
366 fPadPedestalArrayEvent.Delete();
367
368 if ( fDebugStreamer) delete fDebugStreamer;
3cd27a08 369 delete fROC;
467a4b25 370 delete fParam;
371}
372//_____________________________________________________________________
3cd27a08 373Int_t AliTPCCalibPulser::Update(const Int_t icsector,
467a4b25 374 const Int_t icRow,
375 const Int_t icPad,
376 const Int_t icTimeBin,
377 const Float_t csignal)
378{
379 //
3cd27a08 380 // Signal filling methode on the fly pedestal and time offset correction if necessary.
467a4b25 381 // no extra analysis necessary. Assumes knowledge of the signal shape!
382 // assumes that it is looped over consecutive time bins of one pad
383 //
384 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
385
386 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
387
388 //init first pad and sector in this event
389 if ( fCurrentChannel == -1 ) {
390 fCurrentChannel = iChannel;
391 fCurrentSector = icsector;
392 fCurrentRow = icRow;
393 }
394
395 //process last pad if we change to a new one
396 if ( iChannel != fCurrentChannel ){
397 ProcessPad();
398 fCurrentChannel = iChannel;
399 fCurrentSector = icsector;
400 fCurrentRow = icRow;
401 }
402
403 //fill signals for current pad
404 fPadSignal[icTimeBin]=csignal;
405 if ( csignal > fMaxPadSignal ){
406 fMaxPadSignal = csignal;
407 fMaxTimeBin = icTimeBin;
408 }
409 return 0;
410}
411//_____________________________________________________________________
412void AliTPCCalibPulser::FindPedestal(Float_t part)
413{
414 //
415 // find pedestal and noise for the current pad. Use either database or
416 // truncated mean with part*100%
417 //
418 Bool_t noPedestal = kTRUE;;
419 if (fPedestalTPC&&fPadNoiseTPC){
420 //use pedestal database
421 //only load new pedestals if the sector has changed
422 if ( fCurrentSector!=fLastSector ){
423 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
424 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
425 fLastSector=fCurrentSector;
426 }
427
428 if ( fPedestalROC&&fPadNoiseROC ){
429 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
430 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
431 noPedestal = kFALSE;
432 }
433
434 }
435
436 //if we are not running with pedestal database, or for the current sector there is no information
437 //available, calculate the pedestal and noise on the fly
438 if ( noPedestal ) {
439 const Int_t kPedMax = 100; //maximum pedestal value
440 Float_t max = 0;
441 Float_t maxPos = 0;
442 Int_t median = -1;
443 Int_t count0 = 0;
444 Int_t count1 = 0;
445 //
446 Float_t padSignal=0;
447 //
448 UShort_t histo[kPedMax];
449 memset(histo,0,kPedMax*sizeof(UShort_t));
450
451 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
452 padSignal = fPadSignal.GetMatrixArray()[i];
453 if (padSignal<=0) continue;
454 if (padSignal>max && i>10) {
455 max = padSignal;
456 maxPos = i;
457 }
458 if (padSignal>kPedMax-1) continue;
459 histo[Int_t(padSignal+0.5)]++;
460 count0++;
461 }
462 //
463 for (Int_t i=1; i<kPedMax; ++i){
464 if (count1<count0*0.5) median=i;
465 count1+=histo[i];
466 }
467 // truncated mean
468 //
469 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
470 //
471 for (Int_t idelta=1; idelta<10; ++idelta){
472 if (median-idelta<=0) continue;
473 if (median+idelta>kPedMax) continue;
474 if (count<part*count1){
475 count+=histo[median-idelta];
476 mean +=histo[median-idelta]*(median-idelta);
477 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
478 count+=histo[median+idelta];
479 mean +=histo[median+idelta]*(median+idelta);
480 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
481 }
482 }
483 fPadPedestal = 0;
484 fPadNoise = 0;
485 if ( count > 0 ) {
486 mean/=count;
487 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
488 fPadPedestal = mean;
489 fPadNoise = rms;
490 }
491 }
492}
493//_____________________________________________________________________
494void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
495{
496 //
497 // Find position, signal width and height of the CE signal (last signal)
498 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
499 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
500 //
501
502 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
503 Int_t cemaxpos = fMaxTimeBin;
504 Float_t ceSumThreshold = 8.*fPadNoise; // threshold for the signal sum
505 const Int_t kCemin = 2; // range for the analysis of the ce signal +- channels from the peak
506 const Int_t kCemax = 7;
507
508 if (cemaxpos!=0){
509 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
510 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
511 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
512 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
513 ceTime+=signal*(i+0.5);
514 ceRMS +=signal*(i+0.5)*(i+0.5);
515 ceQsum+=signal;
516 }
517 }
518 }
519 if (ceQmax&&ceQsum>ceSumThreshold) {
520 ceTime/=ceQsum;
521 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
522 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
523 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
524
525 //Normalise Q to the pad area
526 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
527
528 ceQsum/=norm;
529 } else {
530 ceQmax=0;
531 ceTime=0;
532 ceRMS =0;
533 ceQsum=0;
534 }
535 param[0] = ceQmax;
536 param[1] = ceTime;
537 param[2] = ceRMS;
538 qSum = ceQsum;
539}
540//_____________________________________________________________________
3cd27a08 541void AliTPCCalibPulser::ProcessPad()
467a4b25 542{
543 //
544 // Process data of current pad
545 //
546
547 FindPedestal();
548 TVectorD param(3);
3cd27a08 549 Float_t qSum;
550 FindPulserSignal(param, qSum);
467a4b25 551
552 Double_t meanT = param[1];
553 Double_t sigmaT = param[2];
554
555 //Fill Event T0 counter
556 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
557
558 //Fill Q histogram
3cd27a08 559 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
467a4b25 560
561 //Fill RMS histogram
562 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
563
564
565 //Fill debugging info
566 if ( fDebugLevel>0 ){
567 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
568 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
3cd27a08 569 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
467a4b25 570 }
467a4b25 571 ResetPad();
572}
573//_____________________________________________________________________
3cd27a08 574void AliTPCCalibPulser::EndEvent()
467a4b25 575{
576 //
577 // Process data of current event
578 //
579
580 //check if last pad has allready been processed, if not do so
581 if ( fMaxTimeBin>-1 ) ProcessPad();
582
583 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
584 for ( Int_t iSec = 0; iSec<72; ++iSec ){
585 TVectorF *vTimes = GetPadTimesEvent(iSec);
586 if ( !vTimes ) continue;
587
588 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
3cd27a08 589 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
590 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
467a4b25 591
3cd27a08 592 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
467a4b25 593
594
595 //Debug start
596 if ( fDebugLevel>0 ){
597 if ( !fDebugStreamer ) {
598 //debug stream
599 TDirectory *backup = gDirectory;
600 fDebugStreamer = new TTreeSRedirector("deb2.root");
601 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
602 }
603
604 Int_t row=0;
605 Int_t pad=0;
606 Int_t padc=0;
607
3cd27a08 608 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
609 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
467a4b25 610
611 UInt_t channel=iChannel;
612 Int_t sector=iSec;
613
614 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
615 pad = channel-fROC->GetRowIndexes(sector)[row];
616 padc = pad-(fROC->GetNPads(sector,row)/2);
617
618 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
619 Form("hSignalD%d.%d.%d",sector,row,pad),
620 fLastTimeBin-fFirstTimeBin,
621 fFirstTimeBin,fLastTimeBin);
622 h1->SetDirectory(0);
623
624 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
625 h1->Fill(i,fPadSignal(i));
626
627 (*fDebugStreamer) << "DataPad" <<
3cd27a08 628// "Event=" << fEvent <<
467a4b25 629 "Sector="<< sector <<
630 "Row=" << row<<
631 "Pad=" << pad <<
632 "PadC=" << padc <<
633 "PadSec="<< channel <<
3cd27a08 634 "Time0=" << time0 <<
635 "Time=" << time <<
636 "RMS=" << rms <<
637 "Sum=" << q <<
467a4b25 638 "hist.=" << h1 <<
639 "\n";
640
641 delete h1;
642 }
643 //Debug end
644
645 }
646 }
647
648}
649//_____________________________________________________________________
3cd27a08 650Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
467a4b25 651{
652 //
653 // Event Processing loop - AliTPCRawStream
654 //
655
656 rawStream->SetOldRCUFormat(fOldRCUformat);
657
658 ResetEvent();
659
660 Bool_t withInput = kFALSE;
661
662 while (rawStream->Next()) {
663
664 Int_t isector = rawStream->GetSector(); // current sector
665 Int_t iRow = rawStream->GetRow(); // current row
666 Int_t iPad = rawStream->GetPad(); // current pad
667 Int_t iTimeBin = rawStream->GetTime(); // current time bin
668 Float_t signal = rawStream->GetSignal(); // current ADC signal
669
670 Update(isector,iRow,iPad,iTimeBin,signal);
671 withInput = kTRUE;
672 }
673
674 if (withInput){
675 EndEvent();
676 }
677
678 return withInput;
679}
680//_____________________________________________________________________
681Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
682{
683 //
684 // Event processing loop - AliRawReader
685 //
686
687
688 AliTPCRawStream rawStream(rawReader);
689
690 rawReader->Select("TPC");
691
692 return ProcessEvent(&rawStream);
693}
694//_____________________________________________________________________
695Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
696{
697 //
698 // Event processing loop - date event
699 //
700 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
701 Bool_t result=ProcessEvent(rawReader);
702 delete rawReader;
703 return result;
704
705}
706//_____________________________________________________________________
3cd27a08 707TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
467a4b25 708 Int_t nbinsY, Float_t ymin, Float_t ymax,
709 Char_t *type, Bool_t force)
710{
711 //
712 // return pointer to Q histogram
713 // if force is true create a new histogram if it doesn't exist allready
714 //
715 if ( !force || arr->UncheckedAt(sector) )
716 return (TH2S*)arr->UncheckedAt(sector);
717
718 // if we are forced and histogram doesn't yes exist create it
719 Char_t name[255], title[255];
720
721 sprintf(name,"hCalib%s%.2d",type,sector);
722 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
723
724 // new histogram with Q calib information. One value for each pad!
725 TH2S* hist = new TH2S(name,title,
726 nbinsY, ymin, ymax,
727 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
728 hist->SetDirectory(0);
729 arr->AddAt(hist,sector);
730 return hist;
731}
732//_____________________________________________________________________
3cd27a08 733TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
467a4b25 734{
735 //
736 // return pointer to T0 histogram
737 // if force is true create a new histogram if it doesn't exist allready
738 //
739 TObjArray *arr = &fHistoT0Array;
740 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
741}
742//_____________________________________________________________________
3cd27a08 743TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
467a4b25 744{
745 //
746 // return pointer to Q histogram
747 // if force is true create a new histogram if it doesn't exist allready
748 //
749 TObjArray *arr = &fHistoQArray;
750 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
751}
752//_____________________________________________________________________
3cd27a08 753TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
467a4b25 754{
755 //
756 // return pointer to Q histogram
757 // if force is true create a new histogram if it doesn't exist allready
758 //
759 TObjArray *arr = &fHistoRMSArray;
760 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
761}
762//_____________________________________________________________________
3cd27a08 763TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
467a4b25 764{
765 //
766 // return pointer to Pad Info from 'arr' for the current event and sector
767 // if force is true create it if it doesn't exist allready
768 //
769 if ( !force || arr->UncheckedAt(sector) )
770 return (TVectorF*)arr->UncheckedAt(sector);
771
772 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
773 arr->AddAt(vect,sector);
774 return vect;
775}
776//_____________________________________________________________________
3cd27a08 777TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
467a4b25 778{
779 //
780 // return pointer to Pad Times Array for the current event and sector
781 // if force is true create it if it doesn't exist allready
782 //
783 TObjArray *arr = &fPadTimesArrayEvent;
784 return GetPadInfoEvent(sector,arr,force);
785}
786//_____________________________________________________________________
3cd27a08 787TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
467a4b25 788{
789 //
790 // return pointer to Pad Q Array for the current event and sector
791 // if force is true create it if it doesn't exist allready
792 // for debugging purposes only
793 //
794
795 TObjArray *arr = &fPadQArrayEvent;
796 return GetPadInfoEvent(sector,arr,force);
797}
798//_____________________________________________________________________
3cd27a08 799TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
467a4b25 800{
801 //
802 // return pointer to Pad RMS Array for the current event and sector
803 // if force is true create it if it doesn't exist allready
804 // for debugging purposes only
805 //
806 TObjArray *arr = &fPadRMSArrayEvent;
807 return GetPadInfoEvent(sector,arr,force);
808}
809//_____________________________________________________________________
3cd27a08 810TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
467a4b25 811{
812 //
813 // return pointer to Pad RMS Array for the current event and sector
814 // if force is true create it if it doesn't exist allready
815 // for debugging purposes only
816 //
817 TObjArray *arr = &fPadPedestalArrayEvent;
818 return GetPadInfoEvent(sector,arr,force);
819}
820//_____________________________________________________________________
3cd27a08 821AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
467a4b25 822{
823 //
824 // return pointer to ROC Calibration
825 // if force is true create a new histogram if it doesn't exist allready
826 //
827 if ( !force || arr->UncheckedAt(sector) )
828 return (AliTPCCalROC*)arr->UncheckedAt(sector);
829
830 // if we are forced and histogram doesn't yes exist create it
831
832 // new AliTPCCalROC for T0 information. One value for each pad!
833 AliTPCCalROC *croc = new AliTPCCalROC(sector);
834 arr->AddAt(croc,sector);
835 return croc;
836}
837//_____________________________________________________________________
3cd27a08 838AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
467a4b25 839{
840 //
841 // return pointer to Carge ROC Calibration
842 // if force is true create a new histogram if it doesn't exist allready
843 //
844 TObjArray *arr = &fCalRocArrayT0;
845 return GetCalRoc(sector, arr, force);
846}
847//_____________________________________________________________________
3cd27a08 848AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
467a4b25 849{
850 //
851 // return pointer to T0 ROC Calibration
852 // if force is true create a new histogram if it doesn't exist allready
853 //
854 TObjArray *arr = &fCalRocArrayQ;
855 return GetCalRoc(sector, arr, force);
856}
857//_____________________________________________________________________
3cd27a08 858AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
467a4b25 859{
860 //
861 // return pointer to signal width ROC Calibration
862 // if force is true create a new histogram if it doesn't exist allready
863 //
864 TObjArray *arr = &fCalRocArrayRMS;
865 return GetCalRoc(sector, arr, force);
866}
867//_____________________________________________________________________
868AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
869{
870 //
871 // return pointer to Outliers
872 // if force is true create a new histogram if it doesn't exist allready
873 //
874 TObjArray *arr = &fCalRocArrayOutliers;
875 return GetCalRoc(sector, arr, force);
876}
877//_____________________________________________________________________
3cd27a08 878void AliTPCCalibPulser::ResetEvent()
467a4b25 879{
880 //
881 // Reset global counters -- Should be called before each event is processed
882 //
883 fLastSector=-1;
884 fCurrentSector=-1;
885 fCurrentRow=-1;
886 fCurrentChannel=-1;
887
888 ResetPad();
889
890 fPadTimesArrayEvent.Delete();
891 fPadQArrayEvent.Delete();
892 fPadRMSArrayEvent.Delete();
893 fPadPedestalArrayEvent.Delete();
894
895 for ( Int_t i=0; i<72; ++i ){
896 fVTime0Offset[i]=0;
897 fVTime0OffsetCounter[i]=0;
898 }
899}
900//_____________________________________________________________________
3cd27a08 901void AliTPCCalibPulser::ResetPad()
467a4b25 902{
903 //
904 // Reset pad infos -- Should be called after a pad has been processed
905 //
906 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
907 fPadSignal[i] = 0;
908 fMaxTimeBin = -1;
909 fMaxPadSignal = -1;
910 fPadPedestal = -1;
911 fPadNoise = -1;
912}
913//_____________________________________________________________________
914void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
915{
916 //
917 // Merge reference histograms of sig to the current AliTPCCalibPulser
918 //
919
920 //merge histograms
921 for (Int_t iSec=0; iSec<72; ++iSec){
922 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
923 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
924 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
925
926
927 if ( hRefQmerge ){
928 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
929 TH2S *hRefQ = GetHistoQ(iSec);
930 if ( hRefQ ) hRefQ->Add(hRefQmerge);
931 else {
932 TH2S *hist = new TH2S(*hRefQmerge);
933 hist->SetDirectory(0);
934 fHistoQArray.AddAt(hist, iSec);
935 }
936 hRefQmerge->SetDirectory(dir);
937 }
938 if ( hRefT0merge ){
939 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
940 TH2S *hRefT0 = GetHistoT0(iSec);
941 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
942 else {
943 TH2S *hist = new TH2S(*hRefT0merge);
944 hist->SetDirectory(0);
945 fHistoT0Array.AddAt(hist, iSec);
946 }
947 hRefT0merge->SetDirectory(dir);
948 }
949 if ( hRefRMSmerge ){
950 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
951 TH2S *hRefRMS = GetHistoRMS(iSec);
952 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
953 else {
954 TH2S *hist = new TH2S(*hRefRMSmerge);
955 hist->SetDirectory(0);
956 fHistoRMSArray.AddAt(hist, iSec);
957 }
958 hRefRMSmerge->SetDirectory(dir);
959 }
960
961 }
962}
963//_____________________________________________________________________
3cd27a08 964void AliTPCCalibPulser::Analyse()
467a4b25 965{
966 //
967 // Calculate calibration constants
968 //
969
970 TVectorD paramQ(3);
971 TVectorD paramT0(3);
972 TVectorD paramRMS(3);
973 TMatrixD dummy(3,3);
974
975 for (Int_t iSec=0; iSec<72; ++iSec){
976 TH2S *hT0 = GetHistoT0(iSec);
977 if (!hT0 ) continue;
978
979 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
980 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
981 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
982 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
983
984 TH2S *hQ = GetHistoQ(iSec);
985 TH2S *hRMS = GetHistoRMS(iSec);
986
3cd27a08 987 Short_t *arrayhQ = hQ->GetArray();
988 Short_t *arrayhT0 = hT0->GetArray();
989 Short_t *arrayhRMS = hRMS->GetArray();
467a4b25 990
991 UInt_t nChannels = fROC->GetNChannels(iSec);
992
993 //debug
994 Int_t row=0;
995 Int_t pad=0;
996 Int_t padc=0;
997 //! debug
998
999 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1000
1001
1002 Float_t cogTime0 = -1000;
1003 Float_t cogQ = -1000;
1004 Float_t cogRMS = -1000;
1005 Float_t cogOut = 0;
1006
1007
1008 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1009 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1010 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1011
1012/*
3cd27a08 1013 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1014 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1015 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
467a4b25 1016 cogQ = paramQ[1];
1017 cogTime0 = paramT0[1];
1018 cogRMS = paramRMS[1];
1019*/
3cd27a08 1020 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1021 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1022 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
467a4b25 1023
1024
1025
1026 /*
1027 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1028 cogOut = 1;
1029 cogTime0 = 0;
1030 cogQ = 0;
1031 cogRMS = 0;
1032 }
1033*/
1034 rocQ->SetValue(iChannel, cogQ*cogQ);
1035 rocT0->SetValue(iChannel, cogTime0);
1036 rocRMS->SetValue(iChannel, cogRMS);
1037 rocOut->SetValue(iChannel, cogOut);
1038
1039
1040 //debug
1041 if ( fDebugLevel > 0 ){
1042 if ( !fDebugStreamer ) {
1043 //debug stream
1044 TDirectory *backup = gDirectory;
1045 fDebugStreamer = new TTreeSRedirector("deb2.root");
1046 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1047 }
1048
1049 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1050 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1051 padc = pad-(fROC->GetNPads(iSec,row)/2);
1052
1053 (*fDebugStreamer) << "DataEnd" <<
1054 "Sector=" << iSec <<
1055 "Pad=" << pad <<
1056 "PadC=" << padc <<
1057 "Row=" << row <<
1058 "PadSec=" << iChannel <<
1059 "Q=" << cogQ <<
1060 "T0=" << cogTime0 <<
1061 "RMS=" << cogRMS <<
1062 "\n";
1063 }
1064 //! debug
1065
1066 }
1067
1068 }
1069 delete fDebugStreamer;
1070 fDebugStreamer = 0x0;
1071}
1072//_____________________________________________________________________
1073void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1074{
1075 //
1076 // Write class to file
1077 //
1078
1079 TString sDir(dir);
1080 TString option;
1081
1082 if ( append )
1083 option = "update";
1084 else
1085 option = "recreate";
1086
1087 TDirectory *backup = gDirectory;
1088 TFile f(filename,option.Data());
1089 f.cd();
1090 if ( !sDir.IsNull() ){
1091 f.mkdir(sDir.Data());
1092 f.cd(sDir);
1093 }
1094 this->Write();
1095 f.Close();
1096
1097 if ( backup ) backup->cd();
1098}
1099//_____________________________________________________________________
1100//_________________________ Test Functions ___________________________
1101//_____________________________________________________________________
1102TObjArray* AliTPCCalibPulser::TestBinning()
1103{
1104 //
1105 // Function to test the binning of the reference histograms
1106 // type: T0, Q or RMS
1107 // mode: 0 - number of filled bins per channel
1108 // 1 - number of empty bins between filled bins in one ROC
1109 // returns TObjArray with the test histograms type*2+mode:
1110 // position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1111
1112
1113 TObjArray *histArray = new TObjArray(6);
1114 const Char_t *type[] = {"T0","Q","RMS"};
1115 Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1116
1117 for (Int_t itype = 0; itype<3; ++itype){
1118 for (Int_t imode=0; imode<2; ++imode){
1119 Int_t icount = itype*2+imode;
1120 histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1121 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1122 72,0,72),
1123 icount);
1124 }
1125 }
1126
1127
1128 TH2S *hRef=0x0;
1129 Short_t *array=0x0;
1130 for (Int_t itype = 0; itype<3; ++itype){
1131 for (Int_t iSec=0; iSec<72; ++iSec){
1132 if ( itype == 0 ) hRef = GetHistoT0(iSec);
1133 if ( itype == 1 ) hRef = GetHistoQ(iSec);
1134 if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1135 if ( hRef == 0x0 ) continue;
1136 array = (hRef->GetArray());
1137 UInt_t nChannels = fROC->GetNChannels(iSec);
1138
1139 Int_t nempty=0;
1140 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1141 Int_t nfilled=0;
1142 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1143 Int_t c1 = 0;
1144 Int_t c2 = 0;
1145 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1146 if ( array[offset+iBin]>0 ) {
1147 nfilled++;
1148 if ( c1 && c2 ) nempty++;
1149 else c1 = 1;
1150 }
1151 else if ( c1 ) c2 = 1;
1152
1153 }
1154 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1155 }
1156 ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1157 }
1158 }
1159 return histArray;
1160}