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