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