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