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