]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCalibCE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
remove obsolete library TPCfast
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibCE.cxx
CommitLineData
75d8233f 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
75d8233f 16/* $Id$ */
17
3cd27a08 18////////////////////////////////////////////////////////////////////////////////////////
19// //
20// Implementation of the TPC Central Electrode calibration //
21// //
22// Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
23// //
24////////////////////////////////////////////////////////////////////////////////////////
7fb602b1 25//
26//
27// *************************************************************************************
28// * Class Description *
29// *************************************************************************************
30//
31/* BEGIN_HTML
32 <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
33 using laser runs.</h4>
34
35 The information retrieved is
36 <ul style="list-style-type: square;">
37 <li>Time arrival from the CE</li>
38 <li>Signal width</li>
39 <li>Signal sum</li>
40 </ul>
41
42<h4>Overview:</h4>
43 <ol style="list-style-type: upper-roman;">
44 <li><a href="#working">Working principle</a></li>
45 <li><a href="#user">User interface for filling data</a></li>
46 <li><a href="#info">Stored information</a></li>
47 </ol>
48
49 <h3><a name="working">I. Working principle</a></h3>
50
51 <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
52 (see below). These in the end call the Update(...) function.</h4>
53
54 <ul style="list-style-type: square;">
55 <li>the Update(...) function:<br />
56 In this function the array fPadSignal is filled with the adc signals between the specified range
57 fFirstTimeBin and fLastTimeBin for the current pad.
58 before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
59 stored in fPadSignal.
60 </li>
61 <ul style="list-style-type: square;">
62 <li>the ProcessPad() function:</li>
63 <ol style="list-style-type: decimal;">
64 <li>Find Pedestal and Noise information</li>
65 <ul style="list-style-type: square;">
66 <li>use database information which has to be set by calling<br />
67 SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
68 <li>if no information from the pedestal data base
69 is available the informaion is calculated on the fly
70 ( see FindPedestal() function )</li>
71 </ul>
72 <li>Find local maxima of the pad signal</li>
73 <ul style="list-style-type: square;">
74 <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
75 have been observed ( see FindLocalMaxima(...) )</li>
76 </ul>
77 <li>Find the CE signal information</li>
78 <ul style="list-style-type: square;">
79 <li>to find the position of the CE signal the Tmean information from the previos event is used
80 as the CE signal the local maximum closest to this Tmean is identified</li>
81 <li>calculate mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
82 the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
83 </ul>
84 <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
85 <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
86 </ol>
87 </ul>
88 </ul>
89
90 <h4>At the end of each event the EndEvent() function is called</h4>
91
92 <ul style="list-style-type: square;">
93 <li>the EndEvent() function:</li>
94 <ul style="list-style-type: square;">
95 <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
96 This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
97 <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
98 <li>calculate Mean Q</li>
99 <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
100 </ul>
101 </ul>
102
103 <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
104 <ul style="list-style-type: square;">
105 <li>the Analyse() function:</li>
106 <ul style="list-style-type: square;">
107 <li>calculate the mean values of T0, RMS, Q for each pad, using
108 the AliMathBase::GetCOG(...) function</li>
109 <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
110 (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
111 </ul>
112 </ul>
113
114 <h3><a name="user">II. User interface for filling data</a></h3>
115
116 <h4>To Fill information one of the following functions can be used:</h4>
117
118 <ul style="list-style-type: none;">
119 <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
120 <ul style="list-style-type: square;">
121 <li>process Date event</li>
122 <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
123 </ul>
124 <br />
125
126 <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
127 <ul style="list-style-type: square;">
128 <li>process AliRawReader event</li>
129 <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
130 </ul>
131 <br />
132
133 <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
134 <ul style="list-style-type: square;">
135 <li>process event from AliTPCRawStream</li>
136 <li>call Update function for signal filling</li>
137 </ul>
138 <br />
139
140 <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
141 iPad, const Int_t iTimeBin, const Float_t signal);</li>
142 <ul style="list-style-type: square;">
143 <li>directly fill signal information (sector, row, pad, time bin, pad)
144 to the reference histograms</li>
145 </ul>
146 </ul>
147
148 <h4>It is also possible to merge two independently taken calibrations using the function</h4>
149
150 <ul style="list-style-type: none;">
151 <li> void Merge(AliTPCCalibSignal *sig)</li>
152 <ul style="list-style-type: square;">
153 <li>copy histograms in 'sig' if they do not exist in this instance</li>
154 <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
155 <li>After merging call Analyse again!</li>
156 </ul>
157 </ul>
158
159
160 <h4>example: filling data using root raw data:</h4>
161 <pre>
162 void fillCE(Char_t *filename)
163 {
164 rawReader = new AliRawReaderRoot(fileName);
165 if ( !rawReader ) return;
166 AliTPCCalibCE *calib = new AliTPCCalibCE;
167 while (rawReader->NextEvent()){
168 calib->ProcessEvent(rawReader);
169 }
170 calib->Analyse();
171 calib->DumpToFile("CEData.root");
172 delete rawReader;
173 delete calib;
174 }
175 </pre>
176
177 <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
178
179 <h4><a name="info:stored">III.1 Stored information</a></h4>
180 <ul style="list-style-type: none;">
181 <li>Histograms:</li>
182 <ul style="list-style-type: none;">
183 <li>For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum])
184 is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
185 stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
186 </ul>
187 <br />
188
189 <li>Calibration Data:</li>
190 <ul style="list-style-type: none;">
191 <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
192 the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
193 fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
194 is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
195 </ul>
196 <br />
197
198 <li>For each event the following information is stored:</li>
199
200 <ul style="list-style-type: square;">
201 <li>event time ( TVectorD fVEventTime )</li>
202 <li>event id ( TVectorD fVEventNumber )</li>
203 <br />
204 <li>mean arrival time for each ROC ( TObjArray fTMeanArrayEvent )</li>
205 <li>mean Q for each ROC ( TObjArray fQMeanArrayEvent )</li>
206 <li>parameters of a plane fit for each ROC ( TObjArray fParamArrayEventPol1 )</li>
207 <li>parameters of a 2D parabola fit for each ROC ( TObjArray fParamArrayEventPol2 )</li>
208 </ul>
209 </ul>
210
211 <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
212 <ul style="list-style-type: none;">
213 <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
214 <ul style="list-style-type: square;">
215 <li>TH2F *GetHistoT0(Int_t sector);</li>
216 <li>TH2F *GetHistoRMS(Int_t sector);</li>
217 <li>TH2F *GetHistoQ(Int_t sector);</li>
218 </ul>
219 <br />
220
221 <li>Accessing the calibration storage objects:</li>
222 <ul style="list-style-type: square;">
223 <li>AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values</li>
224 <li>AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values</li>
225 <li>AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values</li>
226 </ul>
227 <br />
228
229 <li>Accessin the event by event information:</li>
230 <ul style="list-style-type: square;">
231 <li>The event by event information can be displayed using the</li>
232 <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
233 <li>which creates a graph from the specified variables</li>
234 </ul>
235 </ul>
236
237 <h4>example for visualisation:</h4>
238 <pre>
239 //if the file "CEData.root" was created using the above example one could do the following:
240 TFile fileCE("CEData.root")
241 AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
242 ce->GetCalRocT0(0)->Draw("colz");
243 ce->GetCalRocRMS(0)->Draw("colz");
244
245 //or use the AliTPCCalPad functionality:
246 AliTPCCalPad padT0(ped->GetCalPadT0());
247 AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
248 padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
249 padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
250
251 //display event by event information:
252 //Draw mean arrival time as a function of the event time for oroc sector A00
253 ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
254 //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
255 ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");
256 </pre>
257END_HTML */
258//////////////////////////////////////////////////////////////////////////////////////
259
260
3cd27a08 261//Root includes
262#include <TObjArray.h>
78f17711 263#include <TH1.h>
3cd27a08 264#include <TH1F.h>
265#include <TH2S.h>
78f17711 266#include <TF1.h>
3cd27a08 267#include <TString.h>
268#include <TVectorF.h>
269#include <TVectorD.h>
78f17711 270#include <TVector3.h>
3cd27a08 271#include <TMatrixD.h>
272#include <TMath.h>
273#include <TGraph.h>
78f17711 274#include <TGraphErrors.h>
3cd27a08 275#include <TString.h>
ac940b58 276#include <TMap.h>
3cd27a08 277#include <TDirectory.h>
278#include <TSystem.h>
279#include <TFile.h>
7442bceb 280#include <TCollection.h>
78f17711 281#include <TTimeStamp.h>
282#include <TList.h>
283#include <TKey.h>
6e6025f4 284#include <TSpectrum.h>
3cd27a08 285
286//AliRoot includes
4c6d06dc 287#include "AliLog.h"
3cd27a08 288#include "AliRawReader.h"
289#include "AliRawReaderRoot.h"
290#include "AliRawReaderDate.h"
291#include "AliRawEventHeaderBase.h"
292#include "AliTPCRawStream.h"
3cd27a08 293#include "AliTPCCalROC.h"
294#include "AliTPCCalPad.h"
295#include "AliTPCROC.h"
296#include "AliTPCParam.h"
297#include "AliTPCCalibCE.h"
298#include "AliMathBase.h"
78f17711 299#include "AliTPCTransform.h"
300#include "AliTPCLaserTrack.h"
3cd27a08 301#include "TTreeStream.h"
302
78f17711 303#include "AliCDBManager.h"
304#include "AliCDBEntry.h"
3cd27a08 305//date
306#include "event.h"
307ClassImp(AliTPCCalibCE)
308
309
7fb602b1 310AliTPCCalibCE::AliTPCCalibCE() :
880c3382 311 AliTPCCalibRawBase(),
312 fNbinsT0(200),
313 fXminT0(-5),
314 fXmaxT0(5),
315 fNbinsQ(200),
316 fXminQ(1),
317 fXmaxQ(40),
318 fNbinsRMS(100),
319 fXminRMS(0.1),
320 fXmaxRMS(5.1),
321 fPeakDetMinus(2),
322 fPeakDetPlus(3),
323 fPeakIntMinus(2),
324 fPeakIntPlus(2),
325 fNoiseThresholdMax(5.),
326 fNoiseThresholdSum(8.),
327 fIsZeroSuppressed(kFALSE),
328 fLastSector(-1),
329 fSecRejectRatio(.4),
330 fParam(new AliTPCParam),
331 fPedestalTPC(0x0),
332 fPadNoiseTPC(0x0),
333 fPedestalROC(0x0),
334 fPadNoiseROC(0x0),
335 fCalRocArrayT0(72),
336 fCalRocArrayT0Err(72),
337 fCalRocArrayQ(72),
338 fCalRocArrayRMS(72),
339 fCalRocArrayOutliers(72),
340 fHistoQArray(72),
341 fHistoT0Array(72),
342 fHistoRMSArray(72),
343 fMeanT0rms(0),
344 fMeanQrms(0),
345 fMeanRMSrms(0),
346 fHistoTmean(72),
347 fParamArrayEventPol1(72),
348 fParamArrayEventPol2(72),
349 fTMeanArrayEvent(72),
350 fQMeanArrayEvent(72),
351 fVEventTime(1000),
352 fVEventNumber(1000),
353 fVTime0SideA(1000),
354 fVTime0SideC(1000),
880c3382 355 fEventId(-1),
c3066940 356 fOldRunNumber(0),
880c3382 357 fPadTimesArrayEvent(72),
358 fPadQArrayEvent(72),
359 fPadRMSArrayEvent(72),
360 fPadPedestalArrayEvent(72),
361 fCurrentChannel(-1),
362 fCurrentSector(-1),
363 fCurrentRow(-1),
364 fMaxPadSignal(-1),
365 fMaxTimeBin(-1),
2963bcbf 366// fPadSignal(1024),
880c3382 367 fPadPedestal(0),
368 fPadNoise(0),
369 fVTime0Offset(72),
370 fVTime0OffsetCounter(72),
371 fVMeanQ(72),
372 fVMeanQCounter(72),
78f17711 373 fCurrentCETimeRef(0),
374 fProcessOld(kTRUE),
375 fProcessNew(kFALSE),
376 fAnalyseNew(kTRUE),
377 fHnDrift(0x0),
378 fArrHnDrift(100),
379 fTimeBursts(100),
6e6025f4 380 fArrFitGraphs(0x0),
381 fEventInBunch(0)
75d8233f 382{
880c3382 383 //
384 // AliTPCSignal default constructor
385 //
386 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
387 fFirstTimeBin=650;
78f17711 388 fLastTimeBin=1030;
880c3382 389 fParam->Update();
2963bcbf 390 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
6e6025f4 391 for (Int_t i=0;i<14;++i){
78f17711 392 fPeaks[i]=0;
393 fPeakWidths[i]=0;
394 }
395 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
75d8233f 396}
397//_____________________________________________________________________
398AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
880c3382 399 AliTPCCalibRawBase(sig),
400 fNbinsT0(sig.fNbinsT0),
401 fXminT0(sig.fXminT0),
402 fXmaxT0(sig.fXmaxT0),
403 fNbinsQ(sig.fNbinsQ),
404 fXminQ(sig.fXminQ),
405 fXmaxQ(sig.fXmaxQ),
406 fNbinsRMS(sig.fNbinsRMS),
407 fXminRMS(sig.fXminRMS),
408 fXmaxRMS(sig.fXmaxRMS),
409 fPeakDetMinus(sig.fPeakDetMinus),
410 fPeakDetPlus(sig.fPeakDetPlus),
411 fPeakIntMinus(sig.fPeakIntMinus),
412 fPeakIntPlus(sig.fPeakIntPlus),
413 fNoiseThresholdMax(sig.fNoiseThresholdMax),
414 fNoiseThresholdSum(sig.fNoiseThresholdSum),
415 fIsZeroSuppressed(sig.fIsZeroSuppressed),
416 fLastSector(-1),
417 fSecRejectRatio(.4),
418 fParam(new AliTPCParam),
419 fPedestalTPC(0x0),
420 fPadNoiseTPC(0x0),
421 fPedestalROC(0x0),
422 fPadNoiseROC(0x0),
423 fCalRocArrayT0(72),
424 fCalRocArrayT0Err(72),
425 fCalRocArrayQ(72),
426 fCalRocArrayRMS(72),
427 fCalRocArrayOutliers(72),
428 fHistoQArray(72),
429 fHistoT0Array(72),
430 fHistoRMSArray(72),
431 fMeanT0rms(sig.fMeanT0rms),
432 fMeanQrms(sig.fMeanQrms),
433 fMeanRMSrms(sig.fMeanRMSrms),
434 fHistoTmean(72),
435 fParamArrayEventPol1(72),
436 fParamArrayEventPol2(72),
437 fTMeanArrayEvent(72),
438 fQMeanArrayEvent(72),
439 fVEventTime(sig.fVEventTime),
440 fVEventNumber(sig.fVEventNumber),
441 fVTime0SideA(sig.fVTime0SideA),
442 fVTime0SideC(sig.fVTime0SideC),
880c3382 443 fEventId(-1),
c3066940 444 fOldRunNumber(0),
880c3382 445 fPadTimesArrayEvent(72),
446 fPadQArrayEvent(72),
447 fPadRMSArrayEvent(72),
448 fPadPedestalArrayEvent(72),
449 fCurrentChannel(-1),
450 fCurrentSector(-1),
451 fCurrentRow(-1),
452 fMaxPadSignal(-1),
453 fMaxTimeBin(-1),
2963bcbf 454// fPadSignal(1024),
880c3382 455 fPadPedestal(0),
456 fPadNoise(0),
457 fVTime0Offset(72),
458 fVTime0OffsetCounter(72),
459 fVMeanQ(72),
460 fVMeanQCounter(72),
78f17711 461 fCurrentCETimeRef(0),
462 fProcessOld(sig.fProcessOld),
463 fProcessNew(sig.fProcessNew),
464 fAnalyseNew(sig.fAnalyseNew),
465 fHnDrift(0x0),
466 fArrHnDrift(100),
467 fTimeBursts(100),
6e6025f4 468 fArrFitGraphs(0x0),
469 fEventInBunch(0)
75d8233f 470{
ac940b58 471 //
880c3382 472 // AliTPCSignal copy constructor
ac940b58 473 //
2963bcbf 474 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
475
ac940b58 476 for (Int_t iSec = 0; iSec < 72; ++iSec){
477 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
478 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
479 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
480 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
75d8233f 481
ac940b58 482 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
483 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
484 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
75d8233f 485
ac940b58 486 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
487 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
488 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
489 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
75d8233f 490
ac940b58 491 if ( hQ != 0x0 ){
ac940b58 492 TH2S *hNew = new TH2S(*hQ);
493 hNew->SetDirectory(0);
494 fHistoQArray.AddAt(hNew,iSec);
ac940b58 495 }
496 if ( hT0 != 0x0 ){
ac940b58 497 TH2S *hNew = new TH2S(*hT0);
498 hNew->SetDirectory(0);
499 fHistoT0Array.AddAt(hNew,iSec);
ac940b58 500 }
501 if ( hRMS != 0x0 ){
ac940b58 502 TH2S *hNew = new TH2S(*hRMS);
503 hNew->SetDirectory(0);
504 fHistoRMSArray.AddAt(hNew,iSec);
75d8233f 505 }
ac940b58 506 }
75d8233f 507
880c3382 508 //copy fit parameters event by event
ac940b58 509 TObjArray *arr=0x0;
510 for (Int_t iSec=0; iSec<72; ++iSec){
511 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
512 if ( arr ){
513 TObjArray *arrEvents = new TObjArray(arr->GetSize());
514 fParamArrayEventPol1.AddAt(arrEvents, iSec);
515 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
516 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
517 arrEvents->AddAt(new TVectorD(*vec),iEvent);
518 }
7fb602b1 519
ac940b58 520 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
521 if ( arr ){
522 TObjArray *arrEvents = new TObjArray(arr->GetSize());
523 fParamArrayEventPol2.AddAt(arrEvents, iSec);
524 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
525 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
526 arrEvents->AddAt(new TVectorD(*vec),iEvent);
7fb602b1 527 }
528
ac940b58 529 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
530 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
531 if ( vMeanTime )
532 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
533 if ( vMeanQ )
534 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
535 }
536
7fb602b1 537
ac940b58 538 fVEventTime.ResizeTo(sig.fVEventTime);
539 fVEventNumber.ResizeTo(sig.fVEventNumber);
540 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
541 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
7fb602b1 542
ac940b58 543 fParam->Update();
78f17711 544
545 for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
546 TObject *o=sig.fArrHnDrift.UncheckedAt(i);
547 if (o){
548 TObject *newo=o->Clone("fHnDrift");
549 fArrHnDrift.AddAt(newo,i);
550 if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
551 }
552 }
553
554 for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
555 fTimeBursts[i]=sig.fTimeBursts[i];
556 }
557
6e6025f4 558 for (Int_t i=0;i<14;++i){
78f17711 559 fPeaks[i]=sig.fPeaks[i];
560 fPeakWidths[i]=sig.fPeakWidths[i];
561 }
562 if (sig.fArrFitGraphs) {
563 fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
564 fArrFitGraphs->SetOwner();
565 }
566
567 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
568
75d8233f 569}
ac940b58 570//_____________________________________________________________________
571AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
880c3382 572 AliTPCCalibRawBase(),
ac940b58 573 fNbinsT0(200),
574 fXminT0(-5),
575 fXmaxT0(5),
576 fNbinsQ(200),
577 fXminQ(1),
578 fXmaxQ(40),
579 fNbinsRMS(100),
580 fXminRMS(0.1),
581 fXmaxRMS(5.1),
880c3382 582 fPeakDetMinus(2),
583 fPeakDetPlus(3),
584 fPeakIntMinus(2),
585 fPeakIntPlus(2),
ac940b58 586 fNoiseThresholdMax(5.),
587 fNoiseThresholdSum(8.),
588 fIsZeroSuppressed(kFALSE),
589 fLastSector(-1),
590 fSecRejectRatio(.4),
ac940b58 591 fParam(new AliTPCParam),
592 fPedestalTPC(0x0),
593 fPadNoiseTPC(0x0),
594 fPedestalROC(0x0),
595 fPadNoiseROC(0x0),
596 fCalRocArrayT0(72),
597 fCalRocArrayT0Err(72),
598 fCalRocArrayQ(72),
599 fCalRocArrayRMS(72),
600 fCalRocArrayOutliers(72),
601 fHistoQArray(72),
602 fHistoT0Array(72),
603 fHistoRMSArray(72),
604 fMeanT0rms(0),
605 fMeanQrms(0),
606 fMeanRMSrms(0),
607 fHistoTmean(72),
608 fParamArrayEventPol1(72),
609 fParamArrayEventPol2(72),
610 fTMeanArrayEvent(72),
611 fQMeanArrayEvent(72),
880c3382 612 fVEventTime(1000),
613 fVEventNumber(1000),
614 fVTime0SideA(1000),
615 fVTime0SideC(1000),
ac940b58 616 fEventId(-1),
c3066940 617 fOldRunNumber(0),
ac940b58 618 fPadTimesArrayEvent(72),
619 fPadQArrayEvent(72),
620 fPadRMSArrayEvent(72),
621 fPadPedestalArrayEvent(72),
622 fCurrentChannel(-1),
623 fCurrentSector(-1),
624 fCurrentRow(-1),
625 fMaxPadSignal(-1),
626 fMaxTimeBin(-1),
2963bcbf 627// fPadSignal(1024),
ac940b58 628 fPadPedestal(0),
629 fPadNoise(0),
630 fVTime0Offset(72),
631 fVTime0OffsetCounter(72),
632 fVMeanQ(72),
633 fVMeanQCounter(72),
78f17711 634 fCurrentCETimeRef(0),
635 fProcessOld(kTRUE),
636 fProcessNew(kFALSE),
637 fAnalyseNew(kTRUE),
638 fHnDrift(0x0),
639 fArrHnDrift(100),
640 fTimeBursts(100),
6e6025f4 641 fArrFitGraphs(0x0),
642 fEventInBunch(0)
ac940b58 643{
644 //
645 // constructor which uses a tmap as input to set some specific parameters
646 //
880c3382 647 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
648 fFirstTimeBin=650;
78f17711 649 fLastTimeBin=1030;
ac940b58 650 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
651 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
652 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
653 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
654 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
655 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
656 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
657 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
658 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
659 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
660 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
880c3382 661 if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
662 if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
663 if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
664 if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
ac940b58 665 if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
666 if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
667 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
880c3382 668 if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
ac940b58 669 if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
670
78f17711 671 if (config->GetValue("ProcessOld")) fProcessOld = (Bool_t)((TObjString*)config->GetValue("ProcessOld"))->GetString().Atoi();
672 if (config->GetValue("ProcessNew")) fProcessNew = (Bool_t)((TObjString*)config->GetValue("ProcessNew"))->GetString().Atoi();
673 if (config->GetValue("AnalyseNew")) fAnalyseNew = (Bool_t)((TObjString*)config->GetValue("AnalyseNew"))->GetString().Atoi();
674
2963bcbf 675 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
6e6025f4 676 for (Int_t i=0;i<14;++i){
78f17711 677 fPeaks[i]=0;
678 fPeakWidths[i]=0;
679 }
2963bcbf 680
ac940b58 681 fParam->Update();
78f17711 682 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
ac940b58 683}
684
75d8233f 685//_____________________________________________________________________
686AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
687{
688 //
689 // assignment operator
690 //
691 if (&source == this) return *this;
692 new (this) AliTPCCalibCE(source);
693
694 return *this;
695}
696//_____________________________________________________________________
697AliTPCCalibCE::~AliTPCCalibCE()
698{
78f17711 699 //
700 // destructor
701 //
702
703 fCalRocArrayT0.Delete();
704 fCalRocArrayT0Err.Delete();
705 fCalRocArrayQ.Delete();
706 fCalRocArrayRMS.Delete();
707 fCalRocArrayOutliers.Delete();
708
709 fHistoQArray.Delete();
710 fHistoT0Array.Delete();
711 fHistoRMSArray.Delete();
712
713 fHistoTmean.Delete();
714
715 fParamArrayEventPol1.Delete();
716 fParamArrayEventPol2.Delete();
717 fTMeanArrayEvent.Delete();
718 fQMeanArrayEvent.Delete();
719
720 fPadTimesArrayEvent.Delete();
721 fPadQArrayEvent.Delete();
722 fPadRMSArrayEvent.Delete();
723 fPadPedestalArrayEvent.Delete();
724
725 fArrHnDrift.SetOwner();
726 fArrHnDrift.Delete();
727
728 if (fArrFitGraphs){
729 fArrFitGraphs->SetOwner();
730 delete fArrFitGraphs;
731 }
75d8233f 732}
733//_____________________________________________________________________
7fb602b1 734Int_t AliTPCCalibCE::Update(const Int_t icsector,
75d8233f 735 const Int_t icRow,
736 const Int_t icPad,
737 const Int_t icTimeBin,
738 const Float_t csignal)
739{
ac940b58 740 //
741 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
742 // no extra analysis necessary. Assumes knowledge of the signal shape!
743 // assumes that it is looped over consecutive time bins of one pad
744 //
4c6d06dc 745
78f17711 746 if (!fProcessOld) return 0;
ac940b58 747 //temp
4c6d06dc 748
b401648b 749 if (icRow<0) return 0;
750 if (icPad<0) return 0;
751 if (icTimeBin<0) return 0;
ac940b58 752 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
75d8233f 753
ac940b58 754 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
75d8233f 755
ac940b58 756 //init first pad and sector in this event
757 if ( fCurrentChannel == -1 ) {
758 fLastSector=-1;
759 fCurrentChannel = iChannel;
760 fCurrentSector = icsector;
761 fCurrentRow = icRow;
762 }
75d8233f 763
ac940b58 764 //process last pad if we change to a new one
765 if ( iChannel != fCurrentChannel ){
766 ProcessPad();
767 fLastSector=fCurrentSector;
768 fCurrentChannel = iChannel;
769 fCurrentSector = icsector;
770 fCurrentRow = icRow;
771 }
75d8233f 772
ac940b58 773 //fill signals for current pad
2963bcbf 774 fPadSignal[icTimeBin]=csignal;
ac940b58 775 if ( csignal > fMaxPadSignal ){
776 fMaxPadSignal = csignal;
777 fMaxTimeBin = icTimeBin;
778 }
779 return 0;
75d8233f 780}
78f17711 781
782//_____________________________________________________________________
783void AliTPCCalibCE::ProcessBunch(const Int_t sector, const Int_t row, const Int_t pad,
784 const Int_t length, const UInt_t startTimeBin, const UShort_t* signal)
785{
786 //
787 // new filling method to fill the THnSparse histogram
788 //
789
790 //only in new processing mode
791 if (!fProcessNew) return;
6e6025f4 792 //don't use the IROCs and inner part of OROCs
78f17711 793 if (sector<36) return;
6e6025f4 794 if (row<40) return;
78f17711 795 //only bunches with reasonable length
796 if (length<3||length>10) return;
797
798 UShort_t timeBin = (UShort_t)startTimeBin;
799 //skip first laser layer
6e6025f4 800 if (timeBin<280) return;
801
78f17711 802 Double_t timeBurst=SetBurstHnDrift();
803
6e6025f4 804 Int_t cePeak=((sector/18)%2)*7+6;
78f17711 805 //after 1 event setup peak ranges
6e6025f4 806 if (fEventInBunch==1 && fPeaks[cePeak]==0) {
78f17711 807 // set time range
808 fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
809 FindLaserLayers();
810 // set time range
811 fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
6e6025f4 812 fHnDrift->Reset();
78f17711 813 }
814
815 // After the first event only fill every 5th bin in a row with the CE information
816 Int_t padFill=pad;
6e6025f4 817 if (fEventInBunch==0||(fPeaks[cePeak]>100&&TMath::Abs((Short_t)fPeaks[cePeak]-(Short_t)timeBin)<(Short_t)fPeakWidths[cePeak])){
78f17711 818 Int_t mod=5;
819 Int_t n=pad/mod;
820 padFill=mod*n+mod/2;
821 }
822
823 //noise removal
6e6025f4 824 if (!IsPeakInRange(timeBin+length/2,sector)) return;
78f17711 825
826 Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
827 (Double_t)padFill,(Double_t)timeBin,timeBurst};
828
829 for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
830 Float_t sig=(Float_t)signal[iTimeBin];
6e6025f4 831 // if (fPeaks[6]>900&&timeBin>(fPeaks[6]-20)&&sig<20) continue;
832 // if (fPeaks[6]>900&&timeBin<(fPeaks[6]-fPeakWidth[6])&&sig<5) continue;
78f17711 833 x[3]=timeBin;
834 fHnDrift->Fill(x,sig);
835 --timeBin;
836 }
837}
838//_____________________________________________________________________
839void AliTPCCalibCE::FindLaserLayers()
840{
841 //
842 // Find the laser layer positoins
843 //
844
6e6025f4 845 //A-side + C-side
846 for (Int_t iside=0;iside<2;++iside){
847 Int_t add=7*iside;
848 //find CE signal position and width
849 fHnDrift->GetAxis(0)->SetRangeUser(36+iside*18,53+iside*18);
850 TH1D *hproj=fHnDrift->Projection(3);
851 hproj->GetXaxis()->SetRangeUser(700,1030);
852 Int_t maxbin=hproj->GetMaximumBin();
853 Double_t binc=hproj->GetBinCenter(maxbin);
854 hproj->GetXaxis()->SetRangeUser(binc-5,binc+5);
855
856 fPeaks[add+6]=(UShort_t)TMath::Nint(binc);
857 // fPeakWidths[4]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
858 fPeakWidths[add+6]=7;
859
860 hproj->GetXaxis()->SetRangeUser(0,maxbin-10);
861 TSpectrum s(6);
862 s.Search(hproj,2,"goff");
863 Int_t index[6];
864 TMath::Sort(6,s.GetPositionX(),index,kFALSE);
865 for (Int_t i=0; i<6; ++i){
866 fPeaks[i+add]=(UShort_t)TMath::Nint(s.GetPositionX()[index[i]]);
867 fPeakWidths[i+add]=5;
78f17711 868 }
6e6025f4 869
870 //other peaks
871
872// Int_t timepos=fPeaks[4]-2*fPeakWidths[4];
873// Int_t width=100;
874
875// for (Int_t i=3; i>=0; --i){
876// hproj->GetXaxis()->SetRangeUser(timepos-width,timepos);
877// fPeaks[i]=hproj->GetMaximumBin();
878// fPeakWidths[i]=(UShort_t)TMath::Nint(10.);
879// width=250;
880// timepos=fPeaks[i]-width/2;
881// }
882
883// for (Int_t i=add; i<add+7; ++i){
884// printf("Peak: %u +- %u\n",fPeaks[i],fPeakWidths[i]);
885// }
886 //check width and reset peak if >100
887// for (Int_t i=0; i<5; ++i){
888// if (fPeakWidths[i]>100) {
889// fPeaks[i]=0;
890// fPeakWidths[i]=0;
891// }
892// }
893
894 delete hproj;
78f17711 895 }
78f17711 896}
897
75d8233f 898//_____________________________________________________________________
899void AliTPCCalibCE::FindPedestal(Float_t part)
900{
ac940b58 901 //
75d8233f 902 // find pedestal and noise for the current pad. Use either database or
903 // truncated mean with part*100%
ac940b58 904 //
905 Bool_t noPedestal = kTRUE;
7fb602b1 906
907 //use pedestal database if set
ac940b58 908 if (fPedestalTPC&&fPadNoiseTPC){
75d8233f 909 //only load new pedestals if the sector has changed
ac940b58 910 if ( fCurrentSector!=fLastSector ){
911 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
912 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
913 }
75d8233f 914
ac940b58 915 if ( fPedestalROC&&fPadNoiseROC ){
916 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
917 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
918 noPedestal = kFALSE;
75d8233f 919 }
920
ac940b58 921 }
922
75d8233f 923 //if we are not running with pedestal database, or for the current sector there is no information
924 //available, calculate the pedestal and noise on the fly
ac940b58 925 if ( noPedestal ) {
926 fPadPedestal = 0;
927 fPadNoise = 0;
928 if ( fIsZeroSuppressed ) return;
929 const Int_t kPedMax = 100; //maximum pedestal value
930 Float_t max = 0;
931 Float_t maxPos = 0;
932 Int_t median = -1;
933 Int_t count0 = 0;
934 Int_t count1 = 0;
935 //
936 Float_t padSignal=0;
937 //
938 UShort_t histo[kPedMax];
939 memset(histo,0,kPedMax*sizeof(UShort_t));
bf57d87d 940
7fb602b1 941 //fill pedestal histogram
ac940b58 942 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
2963bcbf 943 padSignal = fPadSignal[i];
ac940b58 944 if (padSignal<=0) continue;
945 if (padSignal>max && i>10) {
946 max = padSignal;
947 maxPos = i;
948 }
949 if (padSignal>kPedMax-1) continue;
950 histo[int(padSignal+0.5)]++;
951 count0++;
952 }
7fb602b1 953 //find median
ac940b58 954 for (Int_t i=1; i<kPedMax; ++i){
955 if (count1<count0*0.5) median=i;
956 count1+=histo[i];
957 }
75d8233f 958 // truncated mean
ac940b58 959 //
960 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
961 //
962 for (Int_t idelta=1; idelta<10; ++idelta){
963 if (median-idelta<=0) continue;
964 if (median+idelta>kPedMax) continue;
965 if (count<part*count1){
966 count+=histo[median-idelta];
967 mean +=histo[median-idelta]*(median-idelta);
968 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
969 count+=histo[median+idelta];
970 mean +=histo[median+idelta]*(median+idelta);
971 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
972 }
973 }
974 if ( count > 0 ) {
975 mean/=count;
976 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
977 fPadPedestal = mean;
978 fPadNoise = rms;
75d8233f 979 }
ac940b58 980 }
75d8233f 981}
982//_____________________________________________________________________
ac940b58 983void AliTPCCalibCE::UpdateCETimeRef()
984{
985 // Find the time reference of the last valid CE signal in sector
986 // for irocs of the A-Side the reference of the corresponging OROC is returned
987 // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
988 if ( fLastSector == fCurrentSector ) return;
989 Int_t sector=fCurrentSector;
990 if ( sector < 18 ) sector+=36;
991 fCurrentCETimeRef=0;
992 TVectorF *vtRef = GetTMeanEvents(sector);
993 if ( !vtRef ) return;
994 Int_t vtRefSize= vtRef->GetNrows();
995 if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
996 else vtRefSize=fNevents;
997 while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
998 fCurrentCETimeRef=(*vtRef)[vtRefSize];
999 AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
1000}
1001//_____________________________________________________________________
75d8233f 1002void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
1003{
ac940b58 1004 //
75d8233f 1005 // Find position, signal width and height of the CE signal (last signal)
1006 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
1007 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
ac940b58 1008 //
75d8233f 1009
ac940b58 1010 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
1011 Int_t cemaxpos = 0;
1012 Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise; // threshold for the signal sum
880c3382 1013 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
1014 const Int_t kCemax = fPeakIntPlus;
75d8233f 1015
ac940b58 1016 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
bf57d87d 1017
75d8233f 1018 // find maximum closest to the sector mean from the last event
ac940b58 1019 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
7fb602b1 1020 // get sector mean of last event
ac940b58 1021 Float_t tmean = fCurrentCETimeRef;
1022 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
1023 minDist = tmean-maxima[imax];
1024 cemaxpos = (Int_t)maxima[imax];
75d8233f 1025 }
ac940b58 1026 }
880c3382 1027// printf("L1 phase TB: %f\n",GetL1PhaseTB());
ac940b58 1028 if (cemaxpos!=0){
2963bcbf 1029 ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
880c3382 1030 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
ac940b58 1031 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
2963bcbf 1032 Float_t signal = fPadSignal[i]-fPadPedestal;
ac940b58 1033 if (signal>0) {
1034 ceTime+=signal*(i+0.5);
1035 ceRMS +=signal*(i+0.5)*(i+0.5);
1036 ceQsum+=signal;
1037 }
1038 }
75d8233f 1039 }
ac940b58 1040 }
1041 if (ceQmax&&ceQsum>ceSumThreshold) {
1042 ceTime/=ceQsum;
1043 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
880c3382 1044 ceTime-=GetL1PhaseTB();
ac940b58 1045 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
1046 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
1047
1048 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
1049 // the pick-up signal should scale with the pad area. In addition
1050 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
1051 // ratio 2/3. The pad area we express in cm2. We normalise the signal
1052 // to the OROC signal (factor 2/3 for the IROCs).
1053 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
1054 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
1055
1056 ceQsum/=norm;
1057 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
1058 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
1059 } else {
1060 ceQmax=0;
1061 ceTime=0;
1062 ceRMS =0;
1063 ceQsum=0;
1064 }
1065 param[0] = ceQmax;
880c3382 1066 param[1] = ceTime;
ac940b58 1067 param[2] = ceRMS;
1068 qSum = ceQsum;
75d8233f 1069}
1070//_____________________________________________________________________
3cd27a08 1071Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
75d8233f 1072{
ac940b58 1073 //
2963bcbf 1074 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
1075 // and 'tplus' timebins after 'pos'
ac940b58 1076 //
1077 if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
1078 for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
1079 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
1080 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
1081 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
1082 ++iTime2;
1083 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
1084 }
1085 return kTRUE;
75d8233f 1086}
1087//_____________________________________________________________________
1088void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
1089{
ac940b58 1090 //
75d8233f 1091 // Find local maxima on the pad signal and Histogram them
ac940b58 1092 //
4c6d06dc 1093 Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal
ac940b58 1094 Int_t count = 0;
2963bcbf 1095
1096 for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
1097 if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
1098 if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
ac940b58 1099 if (count<maxima.GetNrows()){
1100 maxima.GetMatrixArray()[count++]=i;
1101 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
2963bcbf 1102 i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin fPeakDetMinus+fPeakDetPlus-1
ac940b58 1103 }
75d8233f 1104 }
ac940b58 1105 }
75d8233f 1106}
1107//_____________________________________________________________________
7fb602b1 1108void AliTPCCalibCE::ProcessPad()
75d8233f 1109{
880c3382 1110 //
1111 // Process data of current pad
1112 //
1113 FindPedestal();
1114
1115 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
7fb602b1 1116 // + central electrode and possibly post peaks from the CE signal
1117 // however if we are on a high noise pad a lot more peaks due to the noise might occur
880c3382 1118 FindLocalMaxima(maxima);
1119 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
1120
1121 UpdateCETimeRef(); // update the time refenrence for the current sector
7442bceb 1122 if ( fCurrentCETimeRef<1e-30 ) return; //return if we don't have time 0 info, eg if only one side has laser
880c3382 1123 TVectorD param(3);
1124 Float_t qSum;
1125 FindCESignal(param, qSum, maxima);
1126
1127 Double_t meanT = param[1];
1128 Double_t sigmaT = param[2];
1129
75d8233f 1130 //Fill Event T0 counter
880c3382 1131 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
1132
75d8233f 1133 //Fill Q histogram
880c3382 1134 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
1135
75d8233f 1136 //Fill RMS histogram
880c3382 1137 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
1138
1139
75d8233f 1140 //Fill debugging info
880c3382 1141 if ( GetStreamLevel()>0 ){
1142 (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
1143 (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
1144 (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
1145 }
1146
1147 ResetPad();
75d8233f 1148}
1149//_____________________________________________________________________
7fb602b1 1150void AliTPCCalibCE::EndEvent()
75d8233f 1151{
ac940b58 1152 // Process data of current pad
1153 // The Functions 'SetTimeStamp' and 'SetRunNumber' should be called
1154 // before the EndEvent function to set the event timestamp and number!!!
1155 // This is automatically done if the ProcessEvent(AliRawReader *rawReader)
1156 // function was called
78f17711 1157 if (!fProcessOld) {
6e6025f4 1158 if (fProcessNew){
1159 ++fNevents;
1160 ++fEventInBunch;
1161 }
78f17711 1162 return;
1163 }
1164
ac940b58 1165 //check if last pad has allready been processed, if not do so
1166 if ( fMaxTimeBin>-1 ) ProcessPad();
75d8233f 1167
ac940b58 1168 AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
4c6d06dc 1169
ac940b58 1170 TVectorD param(3);
1171 TMatrixD dummy(3,3);
7fb602b1 1172// TVectorF vMeanTime(72);
1173// TVectorF vMeanQ(72);
ac940b58 1174 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
1175 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
1176
1177 //find mean time0 offset for side A and C
1178 //use only orocs due to the better statistics
1179 Double_t time0Side[2]; //time0 for side A:0 and C:1
1180 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:1
1181 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
1182 for ( Int_t iSec = 36; iSec<72; ++iSec ){
1183 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
1184 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
1185 }
1186 if ( time0SideCount[0] >0 )
1187 time0Side[0]/=time0SideCount[0];
1188 if ( time0SideCount[1] >0 )
1189 time0Side[1]/=time0SideCount[1];
75d8233f 1190 // end find time0 offset
ac940b58 1191 AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1192 Int_t nSecMeanT=0;
1193 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1194 for ( Int_t iSec = 0; iSec<72; ++iSec ){
1195 AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1196 //find median and then calculate the mean around it
1197 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
1198 if ( !hMeanT ) continue;
1199 //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1200 if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1201 hMeanT->Reset();
1202 AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1203 continue;
1204 }
1205
1206 Double_t entries = hMeanT->GetEffectiveEntries();
1207 Double_t sum = 0;
1208 Short_t *arr = hMeanT->GetArray()+1;
1209 Int_t ibin=0;
1210 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1211 sum+=arr[ibin];
1212 if ( sum>=(entries/2.) ) break;
1213 }
1214 Int_t delta = 4;
1215 Int_t firstBin = fFirstTimeBin+ibin-delta;
1216 Int_t lastBin = fFirstTimeBin+ibin+delta;
1217 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1218 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
1219 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1220
7fb602b1 1221 // check boundaries for ebye info of mean time
ac940b58 1222 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1223 Int_t vSize=vMeanTime->GetNrows();
1224 if ( vSize < fNevents+1 ){
1225 vMeanTime->ResizeTo(vSize+100);
1226 }
880c3382 1227
1228 // store mean time for the readout sides
1229 vSize=fVTime0SideA.GetNrows();
1230 if ( vSize < fNevents+1 ){
1231 fVTime0SideA.ResizeTo(vSize+100);
1232 fVTime0SideC.ResizeTo(vSize+100);
1233 }
1234 fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1235 fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
ac940b58 1236
1237 vMeanTime->GetMatrixArray()[fNevents]=median;
1238 nSecMeanT++;
1239 // end find median
1240
1241 TVectorF *vTimes = GetPadTimesEvent(iSec);
1242 if ( !vTimes ) continue; //continue if no time information for this sector is available
1243
1244 AliTPCCalROC calIrocOutliers(0);
1245 AliTPCCalROC calOrocOutliers(36);
1246
1247 // calculate mean Q of the sector
1248 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1249 vSize=vMeanQ->GetNrows();
1250 if ( vSize < fNevents+1 ){
1251 vMeanQ->ResizeTo(vSize+100);
1252 }
1253 Float_t meanQ = 0;
1254 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1255 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1256
1257 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1258 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
75d8233f 1259
1260 //set values for temporary roc calibration class
ac940b58 1261 if ( iSec < 36 ) {
1262 calIroc->SetValue(iChannel, time);
7442bceb 1263 if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
bf57d87d 1264
ac940b58 1265 } else {
1266 calOroc->SetValue(iChannel, time);
7442bceb 1267 if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
ac940b58 1268 }
75d8233f 1269
ac940b58 1270 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
c3066940 1271 // test that we really found the CE signal reliably
1272 if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1273 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
75d8233f 1274
1275
1276
1277 //------------------------------- Debug start ------------------------------
880c3382 1278 if ( GetStreamLevel()>0 ){
1279 TTreeSRedirector *streamer=GetDebugStreamer();
1280 if (streamer){
1281 Int_t row=0;
1282 Int_t pad=0;
1283 Int_t padc=0;
1284
1285 Float_t q = (*GetPadQEvent(iSec))[iChannel];
1286 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1287
1288 UInt_t channel=iChannel;
1289 Int_t sector=iSec;
1290
1291 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1292 pad = channel-fROC->GetRowIndexes(sector)[row];
1293 padc = pad-(fROC->GetNPads(sector,row)/2);
1294
75d8233f 1295// TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1296// Form("hSignalD%d.%d.%d",sector,row,pad),
1297// fLastTimeBin-fFirstTimeBin,
1298// fFirstTimeBin,fLastTimeBin);
1299// h1->SetDirectory(0);
ac940b58 1300 //
7fb602b1 1301// for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
75d8233f 1302// h1->Fill(i,fPadSignal(i));
880c3382 1303
1304 Double_t t0Sec = 0;
1305 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1306 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1307 Double_t t0Side = time0Side[(iSec/18)%2];
1308 (*streamer) << "DataPad" <<
ac940b58 1309 "Event=" << fNevents <<
1310 "RunNumber=" << fRunNumber <<
1311 "TimeStamp=" << fTimeStamp <<
1312 "Sector="<< sector <<
1313 "Row=" << row<<
1314 "Pad=" << pad <<
1315 "PadC=" << padc <<
1316 "PadSec="<< channel <<
1317 "Time0Sec=" << t0Sec <<
1318 "Time0Side=" << t0Side <<
1319 "Time=" << time <<
1320 "RMS=" << rms <<
1321 "Sum=" << q <<
1322 "MeanQ=" << meanQ <<
880c3382 1323 // "hist.=" << h1 <<
ac940b58 1324 "\n";
880c3382 1325
1326 // delete h1;
1327 }
ac940b58 1328 }
880c3382 1329 //----------------------------- Debug end ------------------------------
ac940b58 1330 }// end channel loop
1331
75d8233f 1332
2963bcbf 1333 //do fitting now only in debug mode
1334 if (GetDebugLevel()>0){
1335 TVectorD paramPol1(3);
1336 TVectorD paramPol2(6);
1337 TMatrixD matPol1(3,3);
1338 TMatrixD matPol2(6,6);
1339 Float_t chi2Pol1=0;
1340 Float_t chi2Pol2=0;
1341
1342 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1343 if ( iSec < 36 ){
1344 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1345 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1346 } else {
1347 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1348 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1349 }
1350
1351 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1352 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1353 }
1354
1355 //------------------------------- Debug start ------------------------------
1356 if ( GetStreamLevel()>0 ){
1357 TTreeSRedirector *streamer=GetDebugStreamer();
1358 if ( streamer ) {
1359 (*streamer) << "DataRoc" <<
1360// "Event=" << fEvent <<
1361 "RunNumber=" << fRunNumber <<
1362 "TimeStamp=" << fTimeStamp <<
1363 "Sector="<< iSec <<
1364 "hMeanT.=" << hMeanT <<
1365 "median=" << median <<
1366 "paramPol1.=" << &paramPol1 <<
1367 "paramPol2.=" << &paramPol2 <<
1368 "matPol1.=" << &matPol1 <<
1369 "matPol2.=" << &matPol2 <<
1370 "chi2Pol1=" << chi2Pol1 <<
1371 "chi2Pol2=" << chi2Pol2 <<
1372 "\n";
1373 }
880c3382 1374 }
ac940b58 1375 }
75d8233f 1376 //------------------------------- Debug end ------------------------------
ac940b58 1377 hMeanT->Reset();
1378 }// end sector loop
4d885988 1379 //return if no sector has a valid mean time
ac940b58 1380 if ( nSecMeanT == 0 ) return;
4d885988 1381
1382
7fb602b1 1383// fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1384// fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
ac940b58 1385 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1386 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1387 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1388 }
1389 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1390 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
75d8233f 1391
78f17711 1392 ++fNevents;
6e6025f4 1393 if (fProcessNew) ++fEventInBunch;
ac940b58 1394 fOldRunNumber = fRunNumber;
75d8233f 1395
ac940b58 1396 delete calIroc;
1397 delete calOroc;
1398 AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
75d8233f 1399}
1400//_____________________________________________________________________
7fb602b1 1401TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
75d8233f 1402 Int_t nbinsY, Float_t ymin, Float_t ymax,
a6e0ebfe 1403 const Char_t *type, Bool_t force)
75d8233f 1404{
1405 //
1406 // return pointer to TH2S histogram of 'type'
1407 // if force is true create a new histogram if it doesn't exist allready
1408 //
1409 if ( !force || arr->UncheckedAt(sector) )
a3b590cf 1410 return (TH2S*)arr->UncheckedAt(sector);
75d8233f 1411
7fb602b1 1412 // if we are forced and histogram doesn't exist yet create it
75d8233f 1413 // new histogram with Q calib information. One value for each pad!
a3b590cf 1414 TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
75d8233f 1415 nbinsY, ymin, ymax,
1416 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1417 hist->SetDirectory(0);
1418 arr->AddAt(hist,sector);
1419 return hist;
1420}
1421//_____________________________________________________________________
7fb602b1 1422TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
75d8233f 1423{
1424 //
1425 // return pointer to T0 histogram
1426 // if force is true create a new histogram if it doesn't exist allready
1427 //
1428 TObjArray *arr = &fHistoT0Array;
1429 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1430}
1431//_____________________________________________________________________
7fb602b1 1432TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
75d8233f 1433{
1434 //
1435 // return pointer to Q histogram
1436 // if force is true create a new histogram if it doesn't exist allready
1437 //
1438 TObjArray *arr = &fHistoQArray;
1439 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1440}
1441//_____________________________________________________________________
7fb602b1 1442TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
75d8233f 1443{
1444 //
1445 // return pointer to Q histogram
1446 // if force is true create a new histogram if it doesn't exist allready
1447 //
1448 TObjArray *arr = &fHistoRMSArray;
1449 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1450}
1451//_____________________________________________________________________
1452TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
a6e0ebfe 1453 const Char_t *type, Bool_t force)
75d8233f 1454{
1455 //
1456 // return pointer to TH1S histogram
1457 // if force is true create a new histogram if it doesn't exist allready
1458 //
1459 if ( !force || arr->UncheckedAt(sector) )
a3b590cf 1460 return (TH1S*)arr->UncheckedAt(sector);
75d8233f 1461
1462 // if we are forced and histogram doesn't yes exist create it
4d885988 1463 // new histogram with calib information. One value for each pad!
a3b590cf 1464 TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
75d8233f 1465 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1466 hist->SetDirectory(0);
1467 arr->AddAt(hist,sector);
1468 return hist;
1469}
1470//_____________________________________________________________________
1471TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1472{
1473 //
1474 // return pointer to Q histogram
1475 // if force is true create a new histogram if it doesn't exist allready
1476 //
1477 TObjArray *arr = &fHistoTmean;
1478 return GetHisto(sector, arr, "LastTmean", force);
1479}
1480//_____________________________________________________________________
3cd27a08 1481TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
75d8233f 1482{
1483 //
1484 // return pointer to Pad Info from 'arr' for the current event and sector
1485 // if force is true create it if it doesn't exist allready
1486 //
1487 if ( !force || arr->UncheckedAt(sector) )
1488 return (TVectorF*)arr->UncheckedAt(sector);
1489
7fb602b1 1490 TVectorF *vect = new TVectorF(size);
75d8233f 1491 arr->AddAt(vect,sector);
1492 return vect;
1493}
1494//_____________________________________________________________________
7fb602b1 1495TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
75d8233f 1496{
1497 //
1498 // return pointer to Pad Times Array for the current event and sector
1499 // if force is true create it if it doesn't exist allready
1500 //
1501 TObjArray *arr = &fPadTimesArrayEvent;
7fb602b1 1502 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1503}
1504//_____________________________________________________________________
7fb602b1 1505TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
75d8233f 1506{
1507 //
1508 // return pointer to Pad Q Array for the current event and sector
1509 // if force is true create it if it doesn't exist allready
1510 // for debugging purposes only
1511 //
1512
1513 TObjArray *arr = &fPadQArrayEvent;
7fb602b1 1514 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1515}
1516//_____________________________________________________________________
7fb602b1 1517TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
75d8233f 1518{
1519 //
1520 // return pointer to Pad RMS Array for the current event and sector
1521 // if force is true create it if it doesn't exist allready
1522 // for debugging purposes only
1523 //
1524 TObjArray *arr = &fPadRMSArrayEvent;
7fb602b1 1525 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1526}
1527//_____________________________________________________________________
7fb602b1 1528TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
75d8233f 1529{
1530 //
1531 // return pointer to Pad RMS Array for the current event and sector
1532 // if force is true create it if it doesn't exist allready
1533 // for debugging purposes only
1534 //
1535 TObjArray *arr = &fPadPedestalArrayEvent;
7fb602b1 1536 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1537}
1538//_____________________________________________________________________
1539TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1540{
1541 //
1542 // return pointer to the EbyE info of the mean arrival time for 'sector'
1543 // if force is true create it if it doesn't exist allready
1544 //
1545 TObjArray *arr = &fTMeanArrayEvent;
1546 return GetVectSector(sector,arr,100,force);
75d8233f 1547}
1548//_____________________________________________________________________
7fb602b1 1549TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1550{
1551 //
1552 // return pointer to the EbyE info of the mean arrival time for 'sector'
1553 // if force is true create it if it doesn't exist allready
1554 //
1555 TObjArray *arr = &fQMeanArrayEvent;
1556 return GetVectSector(sector,arr,100,force);
1557}
1558//_____________________________________________________________________
3cd27a08 1559AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1560{
1561 //
1562 // return pointer to ROC Calibration
1563 // if force is true create a new histogram if it doesn't exist allready
1564 //
1565 if ( !force || arr->UncheckedAt(sector) )
1566 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1567
1568 // if we are forced and histogram doesn't yes exist create it
1569
1570 // new AliTPCCalROC for T0 information. One value for each pad!
1571 AliTPCCalROC *croc = new AliTPCCalROC(sector);
75d8233f 1572 arr->AddAt(croc,sector);
1573 return croc;
1574}
1575//_____________________________________________________________________
7fb602b1 1576AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
75d8233f 1577{
1578 //
ef7f7670 1579 // return pointer to Time 0 ROC Calibration
75d8233f 1580 // if force is true create a new histogram if it doesn't exist allready
1581 //
1582 TObjArray *arr = &fCalRocArrayT0;
1583 return GetCalRoc(sector, arr, force);
1584}
1585//_____________________________________________________________________
ef7f7670 1586AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1587{
1588 //
1589 // return pointer to the error of Time 0 ROC Calibration
1590 // if force is true create a new histogram if it doesn't exist allready
1591 //
1592 TObjArray *arr = &fCalRocArrayT0Err;
1593 return GetCalRoc(sector, arr, force);
1594}
1595//_____________________________________________________________________
7fb602b1 1596AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
75d8233f 1597{
1598 //
1599 // return pointer to T0 ROC Calibration
1600 // if force is true create a new histogram if it doesn't exist allready
1601 //
1602 TObjArray *arr = &fCalRocArrayQ;
1603 return GetCalRoc(sector, arr, force);
1604}
1605//_____________________________________________________________________
7fb602b1 1606AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
75d8233f 1607{
1608 //
1609 // return pointer to signal width ROC Calibration
1610 // if force is true create a new histogram if it doesn't exist allready
1611 //
1612 TObjArray *arr = &fCalRocArrayRMS;
1613 return GetCalRoc(sector, arr, force);
1614}
1615//_____________________________________________________________________
1616AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1617{
1618 //
1619 // return pointer to Outliers
1620 // if force is true create a new histogram if it doesn't exist allready
1621 //
1622 TObjArray *arr = &fCalRocArrayOutliers;
1623 return GetCalRoc(sector, arr, force);
1624}
1625//_____________________________________________________________________
3cd27a08 1626TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1627{
1628 //
1629 // return pointer to TObjArray of fit parameters
1630 // if force is true create a new histogram if it doesn't exist allready
1631 //
1632 if ( !force || arr->UncheckedAt(sector) )
1633 return (TObjArray*)arr->UncheckedAt(sector);
1634
1635 // if we are forced and array doesn't yes exist create it
1636
1637 // new TObjArray for parameters
1638 TObjArray *newArr = new TObjArray;
1639 arr->AddAt(newArr,sector);
1640 return newArr;
1641}
1642//_____________________________________________________________________
1643TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1644{
1645 //
1646 // return pointer to TObjArray of fit parameters from plane fit
1647 // if force is true create a new histogram if it doesn't exist allready
1648 //
1649 TObjArray *arr = &fParamArrayEventPol1;
1650 return GetParamArray(sector, arr, force);
1651}
1652//_____________________________________________________________________
1653TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1654{
1655 //
1656 // return pointer to TObjArray of fit parameters from parabola fit
1657 // if force is true create a new histogram if it doesn't exist allready
1658 //
1659 TObjArray *arr = &fParamArrayEventPol2;
1660 return GetParamArray(sector, arr, force);
1661}
78f17711 1662
1663//_____________________________________________________________________
1664void AliTPCCalibCE::CreateDVhist()
1665{
1666 //
1667 // Setup the THnSparse for the drift velocity determination
1668 //
1669
1670 //HnSparse bins
1671 //roc, row, pad, timebin, timestamp
1672 TTimeStamp begin(2010,01,01,0,0,0);
1673 TTimeStamp end(2035,01,01,0,0,0);
1674 Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1675
1676 Int_t bins[kHnBinsDV] = { 72, 96, 140, 1030, nbinsTime};
1677 Double_t xmin[kHnBinsDV] = { 0., 0., 0., 0., (Double_t)begin.GetSec()};
1678 Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1679
1680 fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1681 fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1682 fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1683 fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1684 fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1685 fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
6e6025f4 1686 fHnDrift->Reset();
78f17711 1687}
1688
75d8233f 1689//_____________________________________________________________________
7fb602b1 1690void AliTPCCalibCE::ResetEvent()
75d8233f 1691{
1692 //
1693 // Reset global counters -- Should be called before each event is processed
1694 //
1695 fLastSector=-1;
1696 fCurrentSector=-1;
1697 fCurrentRow=-1;
1698 fCurrentChannel=-1;
1699
1700 ResetPad();
1701
1702 fPadTimesArrayEvent.Delete();
1703 fPadQArrayEvent.Delete();
1704 fPadRMSArrayEvent.Delete();
1705 fPadPedestalArrayEvent.Delete();
1706
7fb602b1 1707 for ( Int_t i=0; i<72; ++i ){
bf57d87d 1708 fVTime0Offset.GetMatrixArray()[i]=0;
1709 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1710 fVMeanQ.GetMatrixArray()[i]=0;
1711 fVMeanQCounter.GetMatrixArray()[i]=0;
75d8233f 1712 }
1713}
1714//_____________________________________________________________________
7fb602b1 1715void AliTPCCalibCE::ResetPad()
75d8233f 1716{
1717 //
1718 // Reset pad infos -- Should be called after a pad has been processed
1719 //
7fb602b1 1720 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
2963bcbf 1721 fPadSignal[i] = 0;
75d8233f 1722 fMaxTimeBin = -1;
1723 fMaxPadSignal = -1;
1724 fPadPedestal = -1;
1725 fPadNoise = -1;
1726}
1727//_____________________________________________________________________
7442bceb 1728void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
7fb602b1 1729{
7442bceb 1730 //
1731 // Merge ce to the current AliTPCCalibCE
1732 //
78f17711 1733 MergeBase(ce);
1734 Int_t nCEevents = ce->GetNeventsProcessed();
7442bceb 1735
78f17711 1736 if (fProcessOld&&ce->fProcessOld){
7442bceb 1737 //merge histograms
78f17711 1738 for (Int_t iSec=0; iSec<72; ++iSec){
1739 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1740 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1741 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1742
1743
1744 if ( hRefQmerge ){
1745 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1746 TH2S *hRefQ = GetHistoQ(iSec);
1747 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1748 else {
1749 TH2S *hist = new TH2S(*hRefQmerge);
1750 hist->SetDirectory(0);
1751 fHistoQArray.AddAt(hist, iSec);
1752 }
1753 hRefQmerge->SetDirectory(dir);
7442bceb 1754 }
78f17711 1755 if ( hRefT0merge ){
1756 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1757 TH2S *hRefT0 = GetHistoT0(iSec);
1758 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1759 else {
1760 TH2S *hist = new TH2S(*hRefT0merge);
1761 hist->SetDirectory(0);
1762 fHistoT0Array.AddAt(hist, iSec);
1763 }
1764 hRefT0merge->SetDirectory(dir);
7442bceb 1765 }
78f17711 1766 if ( hRefRMSmerge ){
1767 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1768 TH2S *hRefRMS = GetHistoRMS(iSec);
1769 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1770 else {
1771 TH2S *hist = new TH2S(*hRefRMSmerge);
1772 hist->SetDirectory(0);
1773 fHistoRMSArray.AddAt(hist, iSec);
1774 }
1775 hRefRMSmerge->SetDirectory(dir);
7442bceb 1776 }
78f17711 1777
7442bceb 1778 }
1779
7fb602b1 1780 // merge time information
7442bceb 1781
7442bceb 1782
78f17711 1783 for (Int_t iSec=0; iSec<72; ++iSec){
1784 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1785 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1786 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1787 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1788
1789 TObjArray *arrPol1 = 0x0;
1790 TObjArray *arrPol2 = 0x0;
1791 TVectorF *vMeanTime = 0x0;
1792 TVectorF *vMeanQ = 0x0;
1793
7442bceb 1794 //resize arrays
78f17711 1795 if ( arrPol1CE && arrPol2CE ){
1796 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1797 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1798 arrPol1->Expand(fNevents+nCEevents);
1799 arrPol2->Expand(fNevents+nCEevents);
1800 }
1801 if ( vMeanTimeCE && vMeanQCE ){
1802 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1803 vMeanQ = GetQMeanEvents(iSec,kTRUE);
1804 vMeanTime->ResizeTo(fNevents+nCEevents);
1805 vMeanQ->ResizeTo(fNevents+nCEevents);
1806 }
1807
1808 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1809 if ( arrPol1CE && arrPol2CE ){
1810 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1811 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1812 if ( paramPol1 && paramPol2 ){
1813 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1814 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1815 }
1816 }
1817 if ( vMeanTimeCE && vMeanQCE ){
1818 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1819 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1820 }
1821 }
7442bceb 1822 }
1823
1824
78f17711 1825
1826 const TVectorD& eventTimes = ce->fVEventTime;
1827 const TVectorD& eventIds = ce->fVEventNumber;
1828 const TVectorF& time0SideA = ce->fVTime0SideA;
1829 const TVectorF& time0SideC = ce->fVTime0SideC;
1830 fVEventTime.ResizeTo(fNevents+nCEevents);
1831 fVEventNumber.ResizeTo(fNevents+nCEevents);
1832 fVTime0SideA.ResizeTo(fNevents+nCEevents);
1833 fVTime0SideC.ResizeTo(fNevents+nCEevents);
1834
7442bceb 1835 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
78f17711 1836 Double_t evTime = eventTimes.GetMatrixArray()[iEvent];
1837 Double_t evId = eventIds.GetMatrixArray()[iEvent];
1838 Float_t t0SideA = time0SideA.GetMatrixArray()[iEvent];
1839 Float_t t0SideC = time0SideC.GetMatrixArray()[iEvent];
1840
1841 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1842 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1843 fVTime0SideA.GetMatrixArray()[fNevents+iEvent] = t0SideA;
1844 fVTime0SideC.GetMatrixArray()[fNevents+iEvent] = t0SideC;
1845 }
1846 }
1847
1848 if (fProcessNew&&ce->fProcessNew) {
1849 if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1850 AliError("Number of bursts in the instances to merge are different. No merging done!");
1851 } else {
1852 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1853 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1854 THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1855 if (h && hce) h->Add(hce);
1856 else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
7442bceb 1857 }
78f17711 1858 //TODO: What to do with fTimeBursts???
7fb602b1 1859 }
7442bceb 1860 }
1861
7442bceb 1862 fNevents+=nCEevents; //increase event counter
1863}
7fb602b1 1864
7442bceb 1865//_____________________________________________________________________
1866Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1867{
1868 //
1869 // Merge all objects of this type in list
1870 //
7fb602b1 1871
7442bceb 1872 Long64_t nmerged=1;
7fb602b1 1873
7442bceb 1874 TIter next(list);
1875 AliTPCCalibCE *ce=0;
1876 TObject *o=0;
7fb602b1 1877
7442bceb 1878 while ( (o=next()) ){
1879 ce=dynamic_cast<AliTPCCalibCE*>(o);
1880 if (ce){
1881 Merge(ce);
1882 ++nmerged;
7fb602b1 1883 }
7442bceb 1884 }
7fb602b1 1885
7442bceb 1886 return nmerged;
7fb602b1 1887}
7442bceb 1888
7fb602b1 1889//_____________________________________________________________________
1890TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
75d8233f 1891{
ac940b58 1892 //
1893 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
880c3382 1894 // or side (-1: A-Side, -2: C-Side)
ac940b58 1895 // xVariable: 0-event time, 1-event id, 2-internal event counter
1896 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1897 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1898 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1899 // not used for mean time and mean Q )
1900 // for an example see class description at the beginning
1901 //
75d8233f 1902
ac940b58 1903 TVectorD *xVar = 0x0;
1904 TObjArray *aType = 0x0;
1905 Int_t npoints=0;
75d8233f 1906
96bf9029 1907 // sanity checks
1908 if ( (sector<-2) || (sector>71) ) return 0x0; //sector outside valid range
1909 if ( (xVariable<0) || (xVariable>2) ) return 0x0; //invalid x-variable
1910 if ( (fitType<0) || (fitType>3) ) return 0x0; //invalid fit type
1911 if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1912 if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1913 if ( sector<0 && fitType!=2) return 0x0; //for side wise information only mean time is available
880c3382 1914
1915 if (sector>=0){
1916 if ( fitType==0 ){
1917 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1918 aType = &fParamArrayEventPol1;
1919 if ( aType->At(sector)==0x0 ) return 0x0;
1920 }
1921 else if ( fitType==1 ){
1922 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1923 aType = &fParamArrayEventPol2;
1924 if ( aType->At(sector)==0x0 ) return 0x0;
1925 }
75d8233f 1926
880c3382 1927 }
ac940b58 1928 if ( xVariable == 0 ) xVar = &fVEventTime;
1929 if ( xVariable == 1 ) xVar = &fVEventNumber;
1930 if ( xVariable == 2 ) {
1931 xVar = new TVectorD(fNevents);
1932 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1933 }
880c3382 1934
a3b590cf 1935 Double_t *x = new Double_t[fNevents];
1936 Double_t *y = new Double_t[fNevents];
1937
ac940b58 1938 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1939 if ( fitType<2 ){
1940 TObjArray *events = (TObjArray*)(aType->At(sector));
1941 if ( events->GetSize()<=ievent ) break;
1942 TVectorD *v = (TVectorD*)(events->At(ievent));
1943 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1944 } else if (fitType == 2) {
1945 Double_t xValue=(*xVar)[ievent];
880c3382 1946 Double_t yValue=0;
1947 if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1948 else if (sector==-1) yValue=fVTime0SideA(ievent);
1949 else if (sector==-2) yValue=fVTime0SideC(ievent);
ac940b58 1950 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1951 }else if (fitType == 3) {
1952 Double_t xValue=(*xVar)[ievent];
1953 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1954 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
75d8233f 1955 }
ac940b58 1956 }
75d8233f 1957
ac940b58 1958 TGraph *gr = new TGraph(npoints);
75d8233f 1959 //sort xVariable increasing
ac940b58 1960 Int_t *sortIndex = new Int_t[npoints];
a3b590cf 1961 TMath::Sort(npoints,x,sortIndex, kFALSE);
ac940b58 1962 for (Int_t i=0;i<npoints;++i){
1963 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1964 }
75d8233f 1965
1966
ac940b58 1967 if ( xVariable == 2 ) delete xVar;
4ce766eb 1968 delete [] x;
1969 delete [] y;
1970 delete [] sortIndex;
ac940b58 1971 return gr;
75d8233f 1972}
1973//_____________________________________________________________________
1974void AliTPCCalibCE::Analyse()
1975{
880c3382 1976 //
1977 // Calculate calibration constants
1978 //
78f17711 1979
1980 if (fProcessOld){
1981 TVectorD paramQ(3);
1982 TVectorD paramT0(3);
1983 TVectorD paramRMS(3);
1984 TMatrixD dummy(3,3);
1985
1986 Float_t channelCounter=0;
1987 fMeanT0rms=0;
1988 fMeanQrms=0;
1989 fMeanRMSrms=0;
1990
1991 for (Int_t iSec=0; iSec<72; ++iSec){
1992 TH2S *hT0 = GetHistoT0(iSec);
1993 if (!hT0 ) continue;
1994
1995 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1996 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1997 AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1998 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1999 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
2000
2001 TH2S *hQ = GetHistoQ(iSec);
2002 TH2S *hRMS = GetHistoRMS(iSec);
2003
2004 Short_t *arrayhQ = hQ->GetArray();
2005 Short_t *arrayhT0 = hT0->GetArray();
2006 Short_t *arrayhRMS = hRMS->GetArray();
2007
2008 UInt_t nChannels = fROC->GetNChannels(iSec);
2009
2010 //debug
2011 Int_t row=0;
2012 Int_t pad=0;
2013 Int_t padc=0;
2014 //! debug
2015
2016 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
2017
2018
2019 Float_t cogTime0 = -1000;
2020 Float_t cogQ = -1000;
2021 Float_t cogRMS = -1000;
2022 Float_t cogOut = 0;
2023 Float_t rms = 0;
2024 Float_t rmsT0 = 0;
2025
2026
2027 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
2028 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
2029 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
2030
2031 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2032 fMeanQrms+=rms;
2033 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2034 fMeanT0rms+=rmsT0;
2035 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2036 fMeanRMSrms+=rms;
2037 channelCounter++;
2038
2039 /*
2040 //outlier specifications
2041 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2042 cogOut = 1;
2043 cogTime0 = 0;
2044 cogQ = 0;
2045 cogRMS = 0;
2046 }
2047 */
2048 rocQ->SetValue(iChannel, cogQ*cogQ);
2049 rocT0->SetValue(iChannel, cogTime0);
2050 rocT0Err->SetValue(iChannel, rmsT0);
2051 rocRMS->SetValue(iChannel, cogRMS);
2052 rocOut->SetValue(iChannel, cogOut);
2053
2054
2055 //debug
2056 if ( GetStreamLevel() > 0 ){
2057 TTreeSRedirector *streamer=GetDebugStreamer();
2058 if ( streamer ) {
2059
2060 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2061 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2062 padc = pad-(fROC->GetNPads(iSec,row)/2);
2063
2064 (*streamer) << "DataEnd" <<
2065 "Sector=" << iSec <<
2066 "Pad=" << pad <<
2067 "PadC=" << padc <<
2068 "Row=" << row <<
2069 "PadSec=" << iChannel <<
2070 "Q=" << cogQ <<
2071 "T0=" << cogTime0 <<
2072 "RMS=" << cogRMS <<
2073 "\n";
2074 }
2075 }
2076 //! debug
2077
2078 }
2079
2080 }
2081 if ( channelCounter>0 ){
2082 fMeanT0rms/=channelCounter;
2083 fMeanQrms/=channelCounter;
2084 fMeanRMSrms/=channelCounter;
2085 }
2086 // if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2087 // delete fDebugStreamer;
2088 // fDebugStreamer = 0x0;
2089 fVEventTime.ResizeTo(fNevents);
2090 fVEventNumber.ResizeTo(fNevents);
2091 fVTime0SideA.ResizeTo(fNevents);
2092 fVTime0SideC.ResizeTo(fNevents);
2093 }
2094
2095 if (fProcessNew&&fAnalyseNew){
2096 AnalyseTrack();
2097 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2098 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2099 h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2100 }
2101 }
2102}
2103
2104
2105
2106
2107//
2108// New functions that also use the laser tracks
2109//
2110
2111
2112
2113//_____________________________________________________________________
2114void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2115{
2116 //
2117 //Find the local maximums for the projections to each axis
2118 //
880c3382 2119
78f17711 2120 //find laser layer positoins
2121 fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2122 FindLaserLayers();
2123 THnSparse *hProj=fHnDrift;
2124 Double_t posCE[4]={0.,0.,0.,0.};
2125 Double_t widthCE[4]={0.,0.,0.,0.};
880c3382 2126
78f17711 2127// if(fPeaks[4]!=0){
2128 // find central electrode position once more, separately for IROC, OROC, A-, C-Side
880c3382 2129
78f17711 2130 for (Int_t i=0; i<4; ++i){
6e6025f4 2131 Int_t ce=(i/2>0)*7+6;
78f17711 2132 hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2133 TH1 *h=fHnDrift->Projection(3);
6e6025f4 2134 h->GetXaxis()->SetRangeUser(fPeaks[ce]-fPeakWidths[ce],fPeaks[ce]+fPeakWidths[ce]);
78f17711 2135 Int_t nbinMax=h->GetMaximumBin();
2136 Double_t maximum=h->GetMaximum();
2137// Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2138// if (nbinMax<700||maximum<maxExpected) continue;
2139 Double_t xbinMax=h->GetBinCenter(nbinMax);
6e6025f4 2140 TF1 fgaus("gaus","gaus",xbinMax-5,xbinMax+5);
78f17711 2141 fgaus.SetParameters(maximum,xbinMax,2);
2142 fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2143 fgaus.SetParLimits(2,0.2,4.);
2144 h->Fit(&fgaus,"RQN");
2145// Double_t deltaX=4*fgaus.GetParameter(2);
2146// xbinMax=fgaus.GetParameter(1);
2147 delete h;
2148 posCE[i]=fgaus.GetParameter(1);
2149 widthCE[i]=4*fgaus.GetParameter(2);
2150 hProj->GetAxis(0)->SetRangeUser(0,72);
2151 }
2152// }
2153 //Current drift velocity
2154 Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2155// cout<<"vdrift="<<vdrift<<endl;
2156
2157 AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2158 //loop over all entries in the histogram
2159 Int_t coord[5];
2160 for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2161 //get entry position and content
2162 Double_t adc=hProj->GetBinContent(ichunk,coord);
2163
880c3382 2164
78f17711 2165 Int_t sector = coord[0]-1;
2166 Int_t row = coord[1]-1;
2167 Int_t pad = coord[2]-1;
2168 Int_t timeBin= coord[3]-1;
2169 Double_t time = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2170 Int_t side = (sector/18)%2;
2171// return;
2172// fPeaks[4]=(UInt_t)posCE[sector/18];
2173// fPeakWidths[4]=(UInt_t)widthCE[sector/18];
880c3382 2174
78f17711 2175 //cuts
2176 if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2177 if (adc < 5 ) continue;
2178 if (IsEdgePad(sector,row,pad)) continue;
2179// if (!IsPeakInRange(timeBin)) continue;
2180// if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2181// if (isCE&&(sector!=0)) continue;
880c3382 2182
78f17711 2183 Int_t padmin=-2, padmax=2;
2184 Int_t timemin=-2, timemax=2;
2185 Int_t minsumperpad=2;
2186 //CE or laser tracks
2187 Bool_t isCE=kFALSE;
2188 if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2189 isCE=kTRUE;
2190 padmin=0;
2191 padmax=0;
2192 timemin=-3;
2193 timemax=7;
2194 }
880c3382 2195
78f17711 2196 //
2197 // Find local maximum and cogs
2198 //
2199 Bool_t isMaximum=kTRUE;
2200 Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2201 Double_t cogY=0, rmsY=0;
2202 Int_t npart=0;
880c3382 2203
78f17711 2204 // for position calculation use
2205 for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2206 Float_t lxyz[3];
2207 fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2208
2209 for(Int_t itime=timemin;itime<=timemax;++itime){
2210
2211 Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2212 Double_t val=hProj->GetBinContent(a);
2213 totalmass+=val;
2214
2215 tbcm +=(timeBin+itime)*val;
2216 padcm+=(pad+ipad)*val;
2217 cogY +=lxyz[1]*val;
2218
2219 rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2220 rmspad+=(pad+ipad)*(pad+ipad)*val;
2221 rmsY +=lxyz[1]*lxyz[1]*val;
2222
2223 if (val>0) ++npart;
2224 if (val>adc) {
2225 isMaximum=kFALSE;
2226 break;
2227 }
2228 }
2229 if (!isMaximum) break;
2230 }
2231
2232 if (!isMaximum||!npart) continue;
2233 if (totalmass<npart*minsumperpad) continue;
2234 if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2235 //for CE we want only one pad by construction
880c3382 2236
78f17711 2237 tbcm/=totalmass;
2238 padcm/=totalmass;
2239 cogY/=totalmass;
2240
2241 rmstb/=totalmass;
2242 rmspad/=totalmass;
2243 rmsY/=totalmass;
2244
2245 rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2246 rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2247 rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2248
2249 Int_t cog=TMath::Nint(padcm);
2250
2251 // timebin --> z position
2252 Float_t zlength=fParam->GetZLength(side);
2253// Float_t timePos=tbcm+(1000-fPeaks[4])
2254 // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2255 Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2256
2257 // local to global transformation--> x and y positions
2258 Float_t padlxyz[3];
2259 fROC->GetPositionLocal(sector,row,pad,padlxyz);
2260
2261 Double_t gxyz[3]={padlxyz[0],cogY,gz};
2262 Double_t lxyz[3]={padlxyz[0],cogY,gz};
2263 Double_t igxyz[3]={0,0,0};
2264 AliTPCTransform t1;
2265 t1.RotatedGlobal2Global(sector,gxyz);
2266
2267 Double_t mindist=0;
2268 Int_t trackID=-1;
2269 Int_t trackID2=-1;
2270
2271 //find track id and index of the position in the track (row)
2272 Int_t index=0;
2273 if (!isCE){
2274 index=row+(sector>35)*fROC->GetNRows(0);
2275 trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2276 } else {
2277 trackID=336+((sector/18)%2);
2278 index= fROC->GetRowIndexes(sector)[row]+pad; // global pad position in sector
2279 if (sector<36) {
2280 index+=(sector%18)*fROC->GetNChannels(sector);
2281 } else {
2282 index+=18*fROC->GetNChannels(0);
2283 index+=(sector%18)*fROC->GetNChannels(sector);
2284 }
2285 //TODO: find out about the multiple peaks in the CE
2286// mindist=TMath::Abs(fPeaks[4]-tbcm);
2287 mindist=1.;
2288 }
2289
2290 // fill track vectors
2291 if (trackID>0){
2292 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2293 Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2294
880c3382 2295
78f17711 2296// travel time effect of light includes
880c3382 2297
78f17711 2298 Double_t raylength=ltr->GetRayLength();
2299 Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2300 ltr->GetXYZ(globmir);
2301 if(trackID<336){
2302 if(side==0){
2303 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2304 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
a7307087 2305 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
78f17711 2306 }
2307 else {
2308 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2309 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
a7307087 2310 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
78f17711 2311 }
2312 }
880c3382 2313
78f17711 2314 if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2315 ltr->fVecSec->GetMatrixArray()[index]=sector;
2316 ltr->fVecP2->GetMatrixArray()[index]=0;
2317 ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2318 ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2319 ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2320 ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2321 ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2322 ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2323 ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2324// ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2325 }
2326 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2327 ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2328 igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2329 igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2330 igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2331 }
2332
2333
2334 if (fStreamLevel>4){
2335 (*GetDebugStreamer()) << "clusters" <<
2336 "run=" << fRunNumber <<
2337 "timestamp=" << timestamp <<
2338 "burst=" << burst <<
2339 "side=" << side <<
2340 "sec=" << sector <<
2341 "row=" << row <<
2342 "pad=" << pad <<
2343 "padCog=" << cog <<
2344 "timebin=" << timeBin <<
2345 "cogCE=" << posCE[sector/18] <<
2346 "withCE=" << widthCE[sector/18] <<
2347 "index=" << index <<
2348
2349 "padcm=" << padcm <<
2350 "rmspad=" << rmspad <<
2351
2352 "cogtb=" << tbcm <<
2353 "rmstb=" << rmstb <<
2354
2355 "npad=" << npart <<
2356
2357 "lx=" << padlxyz[0]<<
2358 "ly=" << cogY <<
2359 "lypad=" << padlxyz[1]<<
2360 "rmsY=" << rmsY <<
2361
2362 "gx=" << gxyz[0] <<
2363 "gy=" << gxyz[1] <<
2364 "gz=" << gxyz[2] <<
2365
2366 "igx=" << igxyz[0] <<
2367 "igy=" << igxyz[1] <<
2368 "igz=" << igxyz[2] <<
2369
2370 "mind=" << mindist <<
2371 "max=" << adc <<
2372 "trackid=" << trackID <<
2373 "trackid2=" << trackID2 <<
2374 "npart=" << npart <<
2375 "\n";
2376 } // end stream levelmgz.fElements
2377
2378 }
2379
2380}
2381
2382//_____________________________________________________________________
2383void AliTPCCalibCE::AnalyseTrack()
2384{
2385 //
2386 // Analyse the tracks
2387 //
2388
2389
2390 AliTPCLaserTrack::LoadTracks();
2391// AliTPCParam *param=0x0;
2392// //cdb run number
2393// AliCDBManager *man=AliCDBManager::Instance();
2394// if (man->GetDefaultStorage()){
2395// AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2396// if (entry){
2397// entry->SetOwner(kTRUE);
2398// param = (AliTPCParam*)(entry->GetObject()->Clone());
2399// }
2400// }
2401// if (param){
2402// if (fParam) delete fParam;
2403// fParam=param;
2404// } else {
2405// AliError("Could not get updated AliTPCParam from OCDB!!!");
2406// }
2407
2408 //Measured and ideal laser tracks
2409 TObjArray* arrMeasured = SetupMeasured();
2410 TObjArray* arrIdeal = AliTPCLaserTrack::GetTracks();
2411 AddCEtoIdeal(arrIdeal);
2412
2413 //find bursts and loop over them
2414 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2415 Double_t timestamp=fTimeBursts[iburst];
2416 AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2417 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2418 if (!fHnDrift) continue;
2419 UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2420 if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2421 fBinsLastAna[iburst]=entries;
2422
2423 for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2424// if (iburst==0) FindLaserLayers();
2425
2426 //reset laser tracks
2427 ResetMeasured(arrMeasured);
2428
2429 //find clusters and associate to the tracks
2430 FindLocalMaxima(arrMeasured, timestamp, iburst);
2431
2432 //calculate drift velocity
2433 CalculateDV(arrIdeal,arrMeasured,iburst);
2434
2435 //Dump information to file if requested
2436 if (fStreamLevel>2){
6e6025f4 2437 //printf("make tree\n");
78f17711 2438 //laser track information
880c3382 2439
78f17711 2440 for (Int_t itrack=0; itrack<338; ++itrack){
2441 TObject *iltr=arrIdeal->UncheckedAt(itrack);
2442 TObject *mltr=arrMeasured->UncheckedAt(itrack);
2443 (*GetDebugStreamer()) << "tracks" <<
2444 "run=" << fRunNumber <<
2445 "time=" << timestamp <<
2446 "burst="<< iburst <<
2447 "iltr.=" << iltr <<
2448 "mltr.=" << mltr <<
2449 "\n";
2450 }
2451 }
2452 }
2453 if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2454}
2455
2456//_____________________________________________________________________
2457Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2458 const Double_t *peakposloc, Int_t &itrackMin2)
2459{
2460 //
2461 // Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2462 //
2463
2464
2465 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2466 TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2467 TVector3 vDir;
2468 TVector3 vSt;
2469
2470 Int_t firstbeam=0;
2471 Int_t lastbeam=336/2;
2472 if ( (sector/18)%2 ) {
2473 firstbeam=336/2;
2474 lastbeam=336;
2475 }
2476
2477 mindist=1000000;
2478 Int_t itrackMin=-1;
2479 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2480 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2481// if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2482 vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2483 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2484 if (TMath::Abs(deltaZ)>40) continue;
2485 vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2486 vSt.RotateZ(ltr->GetAlpha());
2487 vDir.RotateZ(ltr->GetAlpha());
2488
2489 Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2490
2491 if (dist<mindist){
2492 mindist=dist;
2493 itrackMin=itrack;
2494 }
2495
2496 }
2497 itrackMin2=-1;
2498 Float_t mindist2=10;
2499 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2500 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2501 if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2502
2503 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2504 if (TMath::Abs(deltaZ)>40) continue;
2505
2506 Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2507 if (dist>1) continue;
2508
2509 if (dist<mindist2){
2510 mindist2=dist;
2511 itrackMin2=itrack;
2512 }
2513 }
2514 mindist=mindist2;
2515 return itrackMin2;
2516
2517}
2518
2519//_____________________________________________________________________
2520Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2521{
2522 //
2523 // return true if pad is on the edge of a row
2524 //
2525 Int_t edge1 = 0;
2526 if ( pad == edge1 ) return kTRUE;
2527 Int_t edge2 = fROC->GetNPads(sector,row)-1;
2528 if ( pad == edge2 ) return kTRUE;
2529
2530 return kFALSE;
2531}
2532
2533//_____________________________________________________________________
2534TObjArray* AliTPCCalibCE::SetupMeasured()
2535{
2536 //
2537 // setup array of measured laser tracks and CE
2538 //
2539
2540 TObjArray *arrIdeal = AliTPCLaserTrack::GetTracks();
2541 TObjArray *arrMeasured = new TObjArray(338);
2542 arrMeasured->SetOwner();
2543 for(Int_t itrack=0;itrack<336;itrack++){
2544 arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2545 }
2546
2547 //"tracks" for CE
2548 AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2549 ltrce->SetId(336);
2550 ltrce->SetSide(0);
2551 ltrce->fVecSec=new TVectorD(557568/2);
2552 ltrce->fVecP2=new TVectorD(557568/2);
2553 ltrce->fVecPhi=new TVectorD(557568/2);
2554 ltrce->fVecGX=new TVectorD(557568/2);
2555 ltrce->fVecGY=new TVectorD(557568/2);
2556 ltrce->fVecGZ=new TVectorD(557568/2);
2557 ltrce->fVecLX=new TVectorD(557568/2);
2558 ltrce->fVecLY=new TVectorD(557568/2);
2559 ltrce->fVecLZ=new TVectorD(557568/2);
2560
2561 arrMeasured->AddAt(ltrce,336); //CE A-Side
2562
2563 ltrce=new AliTPCLaserTrack;
2564 ltrce->SetId(337);
2565 ltrce->SetSide(1);
2566 ltrce->fVecSec=new TVectorD(557568/2);
2567 ltrce->fVecP2=new TVectorD(557568/2);
2568 ltrce->fVecPhi=new TVectorD(557568/2);
2569 ltrce->fVecGX=new TVectorD(557568/2);
2570 ltrce->fVecGY=new TVectorD(557568/2);
2571 ltrce->fVecGZ=new TVectorD(557568/2);
2572 ltrce->fVecLX=new TVectorD(557568/2);
2573 ltrce->fVecLY=new TVectorD(557568/2);
2574 ltrce->fVecLZ=new TVectorD(557568/2);
2575 arrMeasured->AddAt(ltrce,337); //CE C-Side
2576
2577 return arrMeasured;
2578}
2579
2580//_____________________________________________________________________
2581void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2582{
2583 //
2584 // reset array of measured laser tracks and CE
2585 //
2586 for(Int_t itrack=0;itrack<338;itrack++){
2587 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2588 ltr->fVecSec->Zero();
2589 ltr->fVecP2->Zero();
2590 ltr->fVecPhi->Zero();
2591 ltr->fVecGX->Zero();
2592 ltr->fVecGY->Zero();
2593 ltr->fVecGZ->Zero();
2594 ltr->fVecLX->Zero();
2595 ltr->fVecLY->Zero();
2596 ltr->fVecLZ->Zero();
2597 }
2598}
2599
2600//_____________________________________________________________________
2601void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2602{
2603 //
2604 // Add ideal CE positions to the ideal track data
2605 //
2606
2607 arr->Expand(338);
2608 //"tracks" for CE
2609 AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2610 ltrceA->SetId(336);
2611 ltrceA->SetSide(0);
2612 ltrceA->fVecSec=new TVectorD(557568/2);
2613 ltrceA->fVecP2=new TVectorD(557568/2);
2614 ltrceA->fVecPhi=new TVectorD(557568/2);
2615 ltrceA->fVecGX=new TVectorD(557568/2);
2616 ltrceA->fVecGY=new TVectorD(557568/2);
2617 ltrceA->fVecGZ=new TVectorD(557568/2);
2618 ltrceA->fVecLX=new TVectorD(557568/2);
2619 ltrceA->fVecLY=new TVectorD(557568/2);
2620 ltrceA->fVecLZ=new TVectorD(557568/2);
2621 arr->AddAt(ltrceA,336); //CE A-Side
2622
2623 AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2624 ltrceC->SetId(337);
2625 ltrceC->SetSide(1);
2626 ltrceC->fVecSec=new TVectorD(557568/2);
2627 ltrceC->fVecP2=new TVectorD(557568/2);
2628 ltrceC->fVecPhi=new TVectorD(557568/2);
2629 ltrceC->fVecGX=new TVectorD(557568/2);
2630 ltrceC->fVecGY=new TVectorD(557568/2);
2631 ltrceC->fVecGZ=new TVectorD(557568/2);
2632 ltrceC->fVecLX=new TVectorD(557568/2);
2633 ltrceC->fVecLY=new TVectorD(557568/2);
2634 ltrceC->fVecLZ=new TVectorD(557568/2);
2635 arr->AddAt(ltrceC,337); //CE C-Side
2636
2637 //Calculate ideal positoins
2638 Float_t gxyz[3];
2639 Float_t lxyz[3];
2640 Int_t channelSideA=0;
2641 Int_t channelSideC=0;
2642 Int_t channelSide=0;
2643 AliTPCLaserTrack *ltrce=0x0;
2644
2645 for (Int_t isector=0; isector<72; ++isector){
2646 Int_t side=((isector/18)%2);
2647 for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2648 for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2649 fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2650 fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2651 if (side==0) {
2652 ltrce=ltrceA;
2653 channelSide=channelSideA;
2654 } else {
2655 ltrce=ltrceC;
2656 channelSide=channelSideC;
2657 }
2658
2659 ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2660 ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2661 ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2662 ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2663 ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2664// ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2665 ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2666 ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2667// ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2668
2669 if (side==0){
2670 ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2671 ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2672 ++channelSideA;
2673 }
2674 else {
2675 ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2676 ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2677 ++channelSideC;
2678 }
2679 }
2680 }
2681 }
2682
2683
2684}
2685
2686//_____________________________________________________________________
2687void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2688{
2689 //
2690 // calculate the drift velocity from the reconstructed clusters associated
2691 // to the ideal laser tracks
2692 // use two different fit scenarios: Separate fit for A- and C-Side
2693 // Common fit for A- and C-Side
2694 //
2695
2696 if (!fArrFitGraphs){
2697 fArrFitGraphs=new TObjArray;
2698 }
2699
2700// static TLinearFitter fdriftA(5,"hyp4");
2701// static TLinearFitter fdriftC(5,"hyp4");
2702// static TLinearFitter fdriftAC(6,"hyp5");
2703 Double_t timestamp=fTimeBursts[burst];
2704
2705 static TLinearFitter fdriftA(4,"hyp3");
2706 static TLinearFitter fdriftC(4,"hyp3");
2707 static TLinearFitter fdriftAC(5,"hyp4");
2708 TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2709
2710 Float_t chi2A = 10;
2711 Float_t chi2C = 10;
2712 Float_t chi2AC = 10;
2713 Int_t npointsA=0;
2714 Int_t npointsC=0;
2715 Int_t npointsAC=0;
2716
2717 Double_t minres[3]={20.,1,0.8};
2718 //----
2719 for(Int_t i=0;i<3;i++){
2720
2721 fdriftA.ClearPoints();
2722 fdriftC.ClearPoints();
2723 fdriftAC.ClearPoints();
2724
2725 chi2A = 10;
2726 chi2C = 10;
2727 chi2AC = 10;
2728 npointsA=0;
2729 npointsC=0;
2730 npointsAC=0;
2731
2732 for (Int_t itrack=0; itrack<338; ++itrack){
2733 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2734 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2735
2736 //-- Exclude the tracks which has the biggest inclanation angle
2737 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2738 Int_t clustercounter=0;
2739 Int_t indexMax=159;
880c3382 2740
78f17711 2741 //-- exclude the low intensity tracks
880c3382 2742
78f17711 2743 for (Int_t index=0; index<indexMax; ++index){
2744
2745 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2746 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2747 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2748
2749 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
880c3382 2750 }
78f17711 2751 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2752 clustercounter=0;
2753
880c3382 2754
78f17711 2755 //-- drift length
2756 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
880c3382 2757
78f17711 2758 if (itrack>335) indexMax=557568/2;
2759 for (Int_t index=0; index<indexMax; ++index){
2760 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2761 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2762 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2763 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
880c3382 2764
78f17711 2765 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2766 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2767 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2768 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
880c3382 2769
78f17711 2770 //cut if no track info available
2771 if (iltr->GetBundle()==0) continue;
2772 if (iR<133||mR<133) continue;
6e6025f4 2773 if(TMath::Abs(mltr->fVecP2->GetMatrixArray()[index])>minres[i]) continue;
78f17711 2774
2775 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2776 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2777
2778 //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2779 Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2780
2781 if (iltr->GetSide()==0){
2782 fdriftA.AddPoint(xxx,mdrift,1);
2783 }else{
2784 fdriftC.AddPoint(xxx,mdrift,1);
880c3382 2785 }
78f17711 2786// Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2787 Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2788 fdriftAC.AddPoint(xxx2,mdrift,1);
2789
2790 }//end index loop
2791 }//end laser track loop
2792
2793 //perform fit
2794 fdriftA.Eval();
2795 fdriftC.Eval();
2796 fdriftAC.Eval();
2797
2798
2799
2800 //get fit values
2801 fdriftA.GetParameters(fitA);
2802 fdriftC.GetParameters(fitC);
2803 fdriftAC.GetParameters(fitAC);
2804
2805 //Parameters: 0 linear offset
2806 // 1 mean drift velocity correction factor
2807 // 2 relative global y gradient
2808 // 3 radial deformation
2809 // 4 IROC/OROC offset
2810
2811// FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2812
2813 for (Int_t itrack=0; itrack<338; ++itrack){
2814 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2815 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2816
2817 //-- Exclude the tracks which has the biggest inclanation angle
2818 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2819 Int_t clustercounter=0;
2820 Int_t indexMax=159;
2821
2822 //-- exclude the low intensity tracks
2823
2824 for (Int_t index=0; index<indexMax; ++index){
2825 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2826 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2827 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2828 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
880c3382 2829 }
78f17711 2830 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2831 clustercounter=0;
2832
2833 //-- drift length
2834 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2835
2836 if (itrack>335) indexMax=557568/2;
2837 for (Int_t index=0; index<indexMax; ++index){
2838 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2839 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2840 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2841 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2842
2843 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2844 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2845 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2846 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2847
2848 //cut if no track info available
2849 if (iR<60||mR<60) continue;
2850
2851 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2852 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2853
2854 TVectorD *v=&fitA;
2855 if (iltr->GetSide()==1) v=&fitC;
2856// Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR)+(*v)[4]*( iltr->fVecSec->GetMatrixArray()[index]>35);
2857 Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2858
2859 mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2860
2861 }
2862 }
2863
2864 fitA.ResizeTo(7);
2865 fitC.ResizeTo(7);
2866 fitAC.ResizeTo(8);
2867
2868//set statistics values
2869
2870 npointsA= fdriftA.GetNpoints();
2871 if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2872 fitA[5]=npointsA;
2873 fitA[6]=chi2A;
2874
2875 npointsC= fdriftC.GetNpoints();
2876 if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2877 fitC[5]=npointsC;
2878 fitC[6]=chi2C;
2879
2880 npointsAC= fdriftAC.GetNpoints();
2881 if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2882 fitAC[5]=npointsAC;
2883 fitAC[6]=chi2AC;
2884
2885 if (fStreamLevel>2){
2886 //laser track information
2887 (*GetDebugStreamer()) << "DriftV" <<
2888 "iter=" << i <<
2889 "run=" << fRunNumber <<
2890 "time=" << timestamp <<
2891 "fitA.=" << &fitA <<
2892 "fitC.=" << &fitC <<
2893 "fitAC.=" << &fitAC <<
2894 "\n";
2895
880c3382 2896
ef7f7670 2897 }
880c3382 2898
2899 }
78f17711 2900//-----
2901
2902
2903 //Parameters: 0 linear offset (global)
2904 // 1 mean drift velocity correction factor
2905 // 2 relative global y gradient
2906 // 3 radial deformation
2907 // 4 IROC/OROC offset
2908 // 5 linear offset relative A-C
2909
2910 //get graphs
2911 TGraphErrors *grA[7];
2912 TGraphErrors *grC[7];
2913 TGraphErrors *grAC[8];
2914 TString names("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_");
2915 TString namesAC("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_;GRAPH_MEAN_DELAYC_LASER_ALL_");
2916
2917 TObjArray *arrNames=names.Tokenize(";");
2918 TObjArray *arrNamesAC=namesAC.Tokenize(";");
2919
2920 //A-Side graphs
2921 for (Int_t i=0; i<7; ++i){
2922 TString grName=arrNames->UncheckedAt(i)->GetName();
2923 grName+="A";
2924 grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2925 if (!grA[i]){
2926 grA[i]=new TGraphErrors;
2927 grA[i]->SetName(grName.Data());
2928 grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2929 fArrFitGraphs->Add(grA[i]);
2930 }
2931// Int_t ipoint=grA[i]->GetN();
2932 Int_t ipoint=burst;
2933 grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2934 grA[i]->SetPointError(ipoint,60,.0001);
2935 if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2936 }
2937
2938 //C-Side graphs
2939 for (Int_t i=0; i<7; ++i){
2940 TString grName=arrNames->UncheckedAt(i)->GetName();
2941 grName+="C";
2942 grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2943 if (!grC[i]){
2944 grC[i]=new TGraphErrors;
2945 grC[i]->SetName(grName.Data());
2946 grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2947 fArrFitGraphs->Add(grC[i]);
2948 }
2949// Int_t ipoint=grC[i]->GetN();
2950 Int_t ipoint=burst;
2951 grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2952 grC[i]->SetPointError(ipoint,60,.0001);
6e6025f4 2953 if (i<4) grC[i]->SetPointError(ipoint,60,fdriftC.GetCovarianceMatrixElement(i,i));
78f17711 2954 }
2955
2956 //AC-Side graphs
2957 for (Int_t i=0; i<8; ++i){
2958 TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2959 grName+="AC";
2960 grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2961 if (!grAC[i]){
2962 grAC[i]=new TGraphErrors;
2963 grAC[i]->SetName(grName.Data());
2964 grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2965 fArrFitGraphs->Add(grAC[i]);
2966 }
2967// Int_t ipoint=grAC[i]->GetN();
2968 Int_t ipoint=burst;
2969 grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2970 grAC[i]->SetPointError(ipoint,60,.0001);
6e6025f4 2971 if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftAC.GetCovarianceMatrixElement(i,i));
78f17711 2972 }
2973
2974 if (fDebugLevel>10){
2975 printf("A side fit parameters:\n");
2976 fitA.Print();
2977 printf("\nC side fit parameters:\n");
2978 fitC.Print();
2979 printf("\nAC side fit parameters:\n");
2980 fitAC.Print();
2981 }
2982 delete arrNames;
2983 delete arrNamesAC;
2984}
2985
2986//_____________________________________________________________________
2987Double_t AliTPCCalibCE::SetBurstHnDrift()
2988{
2989 //
2990 // Create a new THnSparse for the current burst
2991 // return the time of the current burst
2992 //
2993 Int_t i=0;
2994 for(i=0; i<fTimeBursts.GetNrows(); ++i){
2995 if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2996 if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2997 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2998 return fTimeBursts(i);
2999 }
3000 }
3001
3002 CreateDVhist();
3003 fArrHnDrift.AddAt(fHnDrift,i);
3004 fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
6e6025f4 3005 for (i=0;i<14;++i){
3006 fPeaks[i]=0;
3007 fPeakWidths[i]=0;
3008 }
3009 fEventInBunch=0;
78f17711 3010 return fTimeStamp;
3011}
3012
3013//_____________________________________________________________________
3014void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
3015{
3016 //
3017 // Write class to file
3018 // option can be specified in the dir option:
3019 // options:
3020 // name=<objname>: the name of the calibration object in file will be <objname>
3021 // type=<type>: the saving type:
3022 // 0 - write the complte object
3023 // 1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
3024 // 2 - like 2, but in addition delete objects that will most probably not be used for calibration
3025 // 3 - store only calibration output, don't store the reference histograms
3026 // and THnSparse (requires Analyse called before)
3027 //
3028 // NOTE: to read the object back, the ReadFromFile function should be used
3029 //
3030
3031 TString sDir(dir);
3032 TString objName=GetName();
3033 Int_t type=0;
3034
3035 //get options
3036 TObjArray *arr=sDir.Tokenize(",");
3037 TIter next(arr);
3038 TObjString *s;
3039 while ( (s=(TObjString*)next()) ){
3040 TString optString=s->GetString();
3041 optString.Remove(TString::kBoth,' ');
3042 if (optString.BeginsWith("name=")){
3043 objName=optString.ReplaceAll("name=","");
3044 }
3045 if (optString.BeginsWith("type=")){
3046 optString.ReplaceAll("type=","");
3047 type=optString.Atoi();
3048 }
3049 }
09d5920f 3050 delete arr;
78f17711 3051
6e6025f4 3052 if ( type==4 ){
3053 // only for the new algorithm
3054 AliTPCCalibCE ce;
3055 ce.fArrFitGraphs=fArrFitGraphs;
3056 ce.fNevents=fNevents;
3057 ce.fTimeBursts.ResizeTo(fTimeBursts.GetNrows());
3058 ce.fTimeBursts=fTimeBursts;
3059 ce.fProcessNew=kTRUE;
3060 TFile f(filename,"recreate");
3061 ce.Write(objName.Data());
3062 fArrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3063 f.Close();
3064 ce.fArrFitGraphs=0x0;
3065 return;
3066 }
3067
3068
78f17711 3069 if (type==1||type==2) {
3070 //delete most probably not needed stuff
3071 //This requires Analyse to be called after reading the object from file
3072 fCalRocArrayT0.Delete();
3073 fCalRocArrayT0Err.Delete();
3074 fCalRocArrayQ.Delete();
3075 fCalRocArrayRMS.Delete();
3076 fCalRocArrayOutliers.Delete();
3077 }
3078 if (type==2||type==3){
3079 fParamArrayEventPol1.Delete();
3080 fParamArrayEventPol2.Delete();
3081 }
3082
3083 TObjArray histoQArray(72);
3084 TObjArray histoT0Array(72);
3085 TObjArray histoRMSArray(72);
3086 TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3087
3088 //save all large 2D histograms in separte pointers
3089 //to have a smaller memory print when saving the object
3090 if (type==1||type==2||type==3){
3091 for (Int_t i=0; i<72; ++i){
3092 histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3093 histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3094 histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3095 }
3096 fHistoQArray.SetOwner(kFALSE);
3097 fHistoT0Array.SetOwner(kFALSE);
3098 fHistoRMSArray.SetOwner(kFALSE);
3099 fHistoQArray.Clear();
3100 fHistoT0Array.Clear();
3101 fHistoRMSArray.Clear();
3102
3103 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3104 arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3105 }
3106 fArrHnDrift.SetOwner(kFALSE);
3107 fArrHnDrift.Clear();
3108 }
3109
3110
3111 TDirectory *backup = gDirectory;
3112
3113 TFile f(filename,"recreate");
6e6025f4 3114 Write(objName.Data());
78f17711 3115 if (type==1||type==2) {
3116 histoQArray.Write("histoQArray",TObject::kSingleKey);
3117 histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3118 histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3119 arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3120 }
3121
3122 f.Save();
3123 f.Close();
3124
3125 //move histograms back to the object
3126 if (type==1||type==2){
3127 for (Int_t i=0; i<72; ++i){
3128 fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3129 fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3130 fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3131 }
3132 fHistoQArray.SetOwner(kTRUE);
3133 fHistoT0Array.SetOwner(kTRUE);
3134 fHistoRMSArray.SetOwner(kTRUE);
3135
3136 for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3137 fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3138 }
3139 fArrHnDrift.SetOwner(kTRUE);
3140 }
3141
3142 if ( backup ) backup->cd();
3143}
3144//_____________________________________________________________________
3145AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3146{
3147 //
3148 // Read object from file
3149 // Handle properly if the histogram arrays were stored separately
3150 // call Analyse to make sure to have the calibration relevant information in the object
3151 //
3152
3153 TFile f(filename);
3154 if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3155 TList *l=f.GetListOfKeys();
3156 TIter next(l);
3157 TKey *key=0x0;
3158 TObject *o=0x0;
3159
3160 AliTPCCalibCE *ce=0x0;
3161 TObjArray *histoQArray=0x0;
3162 TObjArray *histoT0Array=0x0;
3163 TObjArray *histoRMSArray=0x0;
3164 TObjArray *arrHnDrift=0x0;
3165
3166 while ( (key=(TKey*)next()) ){
3167 o=key->ReadObj();
3168 if ( o->IsA()==AliTPCCalibCE::Class() ){
3169 ce=(AliTPCCalibCE*)o;
3170 } else if ( o->IsA()==TObjArray::Class() ){
3171 TString name=key->GetName();
3172 if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3173 if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3174 if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3175 if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3176 }
3177 }
3178
3179 if (ce){
3180 //move histograms back to the object
3181 TH1* hist=0x0;
3182 if (histoQArray){
3183 for (Int_t i=0; i<72; ++i){
3184 hist=(TH1*)histoQArray->UncheckedAt(i);
3185 if (hist) hist->SetDirectory(0x0);
3186 ce->fHistoQArray.AddAt(hist,i);
3187 }
3188 ce->fHistoQArray.SetOwner(kTRUE);
3189 }
3190
3191 if (histoT0Array) {
3192 for (Int_t i=0; i<72; ++i){
3193 hist=(TH1*)histoT0Array->UncheckedAt(i);
3194 if (hist) hist->SetDirectory(0x0);
3195 ce->fHistoT0Array.AddAt(hist,i);
3196 }
3197 ce->fHistoT0Array.SetOwner(kTRUE);
3198 }
3199
3200 if (histoRMSArray){
3201 for (Int_t i=0; i<72; ++i){
3202 hist=(TH1*)histoRMSArray->UncheckedAt(i);
3203 if (hist) hist->SetDirectory(0x0);
3204 ce->fHistoRMSArray.AddAt(hist,i);
3205 }
3206 ce->fHistoRMSArray.SetOwner(kTRUE);
3207 }
3208
3209 if (arrHnDrift){
3210 for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3211 THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3212 ce->fArrHnDrift.AddAt(hSparse,i);
3213 }
3214 }
3215
3216 ce->Analyse();
880c3382 3217 }
78f17711 3218 f.Close();
3219
3220 return ce;
75d8233f 3221}