1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////////////////
20 // Implementation of the TPC Central Electrode calibration //
22 // Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
24 ////////////////////////////////////////////////////////////////////////////////////////
27 // *************************************************************************************
28 // * Class Description *
29 // *************************************************************************************
32 <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
33 using laser runs.</h4>
35 The information retrieved is
36 <ul style="list-style-type: square;">
37 <li>Time arrival from the CE</li>
43 <ol style="list-style-type: upper-roman;">
44 <li><a href="#working">Working principle</a></li>
45 <li><a href="#user">User interface for filling data</a></li>
46 <li><a href="#info">Stored information</a></li>
49 <h3><a name="working">I. Working principle</a></h3>
51 <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
52 (see below). These in the end call the Update(...) function.</h4>
54 <ul style="list-style-type: square;">
55 <li>the Update(...) function:<br />
56 In this function the array fPadSignal is filled with the adc signals between the specified range
57 fFirstTimeBin and fLastTimeBin for the current pad.
58 before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
61 <ul style="list-style-type: square;">
62 <li>the ProcessPad() function:</li>
63 <ol style="list-style-type: decimal;">
64 <li>Find Pedestal and Noise information</li>
65 <ul style="list-style-type: square;">
66 <li>use database information which has to be set by calling<br />
67 SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
68 <li>if no information from the pedestal data base
69 is available the informaion is calculated on the fly
70 ( see FindPedestal() function )</li>
72 <li>Find local maxima of the pad signal</li>
73 <ul style="list-style-type: square;">
74 <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
75 have been observed ( see FindLocalMaxima(...) )</li>
77 <li>Find the CE signal information</li>
78 <ul style="list-style-type: square;">
79 <li>to find the position of the CE signal the Tmean information from the previos event is used
80 as the CE signal the local maximum closest to this Tmean is identified</li>
81 <li>calculate mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
82 the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
84 <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
85 <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
90 <h4>At the end of each event the EndEvent() function is called</h4>
92 <ul style="list-style-type: square;">
93 <li>the EndEvent() function:</li>
94 <ul style="list-style-type: square;">
95 <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
96 This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
97 <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
98 <li>calculate Mean Q</li>
99 <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
103 <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
104 <ul style="list-style-type: square;">
105 <li>the Analyse() function:</li>
106 <ul style="list-style-type: square;">
107 <li>calculate the mean values of T0, RMS, Q for each pad, using
108 the AliMathBase::GetCOG(...) function</li>
109 <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
110 (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
114 <h3><a name="user">II. User interface for filling data</a></h3>
116 <h4>To Fill information one of the following functions can be used:</h4>
118 <ul style="list-style-type: none;">
119 <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
120 <ul style="list-style-type: square;">
121 <li>process Date event</li>
122 <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
126 <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
127 <ul style="list-style-type: square;">
128 <li>process AliRawReader event</li>
129 <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
133 <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
134 <ul style="list-style-type: square;">
135 <li>process event from AliTPCRawStream</li>
136 <li>call Update function for signal filling</li>
140 <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
141 iPad, const Int_t iTimeBin, const Float_t signal);</li>
142 <ul style="list-style-type: square;">
143 <li>directly fill signal information (sector, row, pad, time bin, pad)
144 to the reference histograms</li>
148 <h4>It is also possible to merge two independently taken calibrations using the function</h4>
150 <ul style="list-style-type: none;">
151 <li> void Merge(AliTPCCalibSignal *sig)</li>
152 <ul style="list-style-type: square;">
153 <li>copy histograms in 'sig' if they do not exist in this instance</li>
154 <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
155 <li>After merging call Analyse again!</li>
160 <h4>example: filling data using root raw data:</h4>
162 void fillCE(Char_t *filename)
164 rawReader = new AliRawReaderRoot(fileName);
165 if ( !rawReader ) return;
166 AliTPCCalibCE *calib = new AliTPCCalibCE;
167 while (rawReader->NextEvent()){
168 calib->ProcessEvent(rawReader);
171 calib->DumpToFile("CEData.root");
177 <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
179 <h4><a name="info:stored">III.1 Stored information</a></h4>
180 <ul style="list-style-type: none;">
182 <ul style="list-style-type: none;">
183 <li>For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum])
184 is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
185 stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
189 <li>Calibration Data:</li>
190 <ul style="list-style-type: none;">
191 <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
192 the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
193 fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
194 is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
198 <li>For each event the following information is stored:</li>
200 <ul style="list-style-type: square;">
201 <li>event time ( TVectorD fVEventTime )</li>
202 <li>event id ( TVectorD fVEventNumber )</li>
204 <li>mean arrival time for each ROC ( TObjArray fTMeanArrayEvent )</li>
205 <li>mean Q for each ROC ( TObjArray fQMeanArrayEvent )</li>
206 <li>parameters of a plane fit for each ROC ( TObjArray fParamArrayEventPol1 )</li>
207 <li>parameters of a 2D parabola fit for each ROC ( TObjArray fParamArrayEventPol2 )</li>
211 <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
212 <ul style="list-style-type: none;">
213 <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
214 <ul style="list-style-type: square;">
215 <li>TH2F *GetHistoT0(Int_t sector);</li>
216 <li>TH2F *GetHistoRMS(Int_t sector);</li>
217 <li>TH2F *GetHistoQ(Int_t sector);</li>
221 <li>Accessing the calibration storage objects:</li>
222 <ul style="list-style-type: square;">
223 <li>AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values</li>
224 <li>AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values</li>
225 <li>AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values</li>
229 <li>Accessin the event by event information:</li>
230 <ul style="list-style-type: square;">
231 <li>The event by event information can be displayed using the</li>
232 <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
233 <li>which creates a graph from the specified variables</li>
237 <h4>example for visualisation:</h4>
239 //if the file "CEData.root" was created using the above example one could do the following:
240 TFile fileCE("CEData.root")
241 AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
242 ce->GetCalRocT0(0)->Draw("colz");
243 ce->GetCalRocRMS(0)->Draw("colz");
245 //or use the AliTPCCalPad functionality:
246 AliTPCCalPad padT0(ped->GetCalPadT0());
247 AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
248 padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
249 padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
251 //display event by event information:
252 //Draw mean arrival time as a function of the event time for oroc sector A00
253 ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
254 //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
255 ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");
258 //////////////////////////////////////////////////////////////////////////////////////
262 #include <TObjArray.h>
268 #include <TVectorF.h>
269 #include <TVectorD.h>
270 #include <TVector3.h>
271 #include <TMatrixD.h>
274 #include <TGraphErrors.h>
277 #include <TDirectory.h>
280 #include <TCollection.h>
281 #include <TTimeStamp.h>
287 #include "AliRawReader.h"
288 #include "AliRawReaderRoot.h"
289 #include "AliRawReaderDate.h"
290 #include "AliRawEventHeaderBase.h"
291 #include "AliTPCRawStream.h"
292 #include "AliTPCCalROC.h"
293 #include "AliTPCCalPad.h"
294 #include "AliTPCROC.h"
295 #include "AliTPCParam.h"
296 #include "AliTPCCalibCE.h"
297 #include "AliMathBase.h"
298 #include "AliTPCTransform.h"
299 #include "AliTPCLaserTrack.h"
300 #include "TTreeStream.h"
302 #include "AliCDBManager.h"
303 #include "AliCDBEntry.h"
306 ClassImp(AliTPCCalibCE)
309 AliTPCCalibCE::AliTPCCalibCE() :
310 AliTPCCalibRawBase(),
324 fNoiseThresholdMax(5.),
325 fNoiseThresholdSum(8.),
326 fIsZeroSuppressed(kFALSE),
329 fParam(new AliTPCParam),
335 fCalRocArrayT0Err(72),
338 fCalRocArrayOutliers(72),
346 fParamArrayEventPol1(72),
347 fParamArrayEventPol2(72),
348 fTMeanArrayEvent(72),
349 fQMeanArrayEvent(72),
356 fPadTimesArrayEvent(72),
358 fPadRMSArrayEvent(72),
359 fPadPedestalArrayEvent(72),
369 fVTime0OffsetCounter(72),
372 fCurrentCETimeRef(0),
382 // AliTPCSignal default constructor
384 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
388 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
389 for (Int_t i=0;i<5;++i){
393 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
395 //_____________________________________________________________________
396 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
397 AliTPCCalibRawBase(sig),
398 fNbinsT0(sig.fNbinsT0),
399 fXminT0(sig.fXminT0),
400 fXmaxT0(sig.fXmaxT0),
401 fNbinsQ(sig.fNbinsQ),
404 fNbinsRMS(sig.fNbinsRMS),
405 fXminRMS(sig.fXminRMS),
406 fXmaxRMS(sig.fXmaxRMS),
407 fPeakDetMinus(sig.fPeakDetMinus),
408 fPeakDetPlus(sig.fPeakDetPlus),
409 fPeakIntMinus(sig.fPeakIntMinus),
410 fPeakIntPlus(sig.fPeakIntPlus),
411 fNoiseThresholdMax(sig.fNoiseThresholdMax),
412 fNoiseThresholdSum(sig.fNoiseThresholdSum),
413 fIsZeroSuppressed(sig.fIsZeroSuppressed),
416 fParam(new AliTPCParam),
422 fCalRocArrayT0Err(72),
425 fCalRocArrayOutliers(72),
429 fMeanT0rms(sig.fMeanT0rms),
430 fMeanQrms(sig.fMeanQrms),
431 fMeanRMSrms(sig.fMeanRMSrms),
433 fParamArrayEventPol1(72),
434 fParamArrayEventPol2(72),
435 fTMeanArrayEvent(72),
436 fQMeanArrayEvent(72),
437 fVEventTime(sig.fVEventTime),
438 fVEventNumber(sig.fVEventNumber),
439 fVTime0SideA(sig.fVTime0SideA),
440 fVTime0SideC(sig.fVTime0SideC),
443 fPadTimesArrayEvent(72),
445 fPadRMSArrayEvent(72),
446 fPadPedestalArrayEvent(72),
456 fVTime0OffsetCounter(72),
459 fCurrentCETimeRef(0),
460 fProcessOld(sig.fProcessOld),
461 fProcessNew(sig.fProcessNew),
462 fAnalyseNew(sig.fAnalyseNew),
469 // AliTPCSignal copy constructor
471 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
473 for (Int_t iSec = 0; iSec < 72; ++iSec){
474 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
475 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
476 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
477 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
479 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
480 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
481 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
483 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
484 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
485 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
486 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
489 TH2S *hNew = new TH2S(*hQ);
490 hNew->SetDirectory(0);
491 fHistoQArray.AddAt(hNew,iSec);
494 TH2S *hNew = new TH2S(*hT0);
495 hNew->SetDirectory(0);
496 fHistoT0Array.AddAt(hNew,iSec);
499 TH2S *hNew = new TH2S(*hRMS);
500 hNew->SetDirectory(0);
501 fHistoRMSArray.AddAt(hNew,iSec);
505 //copy fit parameters event by event
507 for (Int_t iSec=0; iSec<72; ++iSec){
508 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
510 TObjArray *arrEvents = new TObjArray(arr->GetSize());
511 fParamArrayEventPol1.AddAt(arrEvents, iSec);
512 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
513 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
514 arrEvents->AddAt(new TVectorD(*vec),iEvent);
517 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
519 TObjArray *arrEvents = new TObjArray(arr->GetSize());
520 fParamArrayEventPol2.AddAt(arrEvents, iSec);
521 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
522 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
523 arrEvents->AddAt(new TVectorD(*vec),iEvent);
526 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
527 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
529 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
531 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
535 fVEventTime.ResizeTo(sig.fVEventTime);
536 fVEventNumber.ResizeTo(sig.fVEventNumber);
537 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
538 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
542 for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
543 TObject *o=sig.fArrHnDrift.UncheckedAt(i);
545 TObject *newo=o->Clone("fHnDrift");
546 fArrHnDrift.AddAt(newo,i);
547 if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
551 for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
552 fTimeBursts[i]=sig.fTimeBursts[i];
555 for (Int_t i=0;i<5;++i){
556 fPeaks[i]=sig.fPeaks[i];
557 fPeakWidths[i]=sig.fPeakWidths[i];
559 if (sig.fArrFitGraphs) {
560 fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
561 fArrFitGraphs->SetOwner();
564 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
567 //_____________________________________________________________________
568 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
569 AliTPCCalibRawBase(),
583 fNoiseThresholdMax(5.),
584 fNoiseThresholdSum(8.),
585 fIsZeroSuppressed(kFALSE),
588 fParam(new AliTPCParam),
594 fCalRocArrayT0Err(72),
597 fCalRocArrayOutliers(72),
605 fParamArrayEventPol1(72),
606 fParamArrayEventPol2(72),
607 fTMeanArrayEvent(72),
608 fQMeanArrayEvent(72),
615 fPadTimesArrayEvent(72),
617 fPadRMSArrayEvent(72),
618 fPadPedestalArrayEvent(72),
628 fVTime0OffsetCounter(72),
631 fCurrentCETimeRef(0),
641 // constructor which uses a tmap as input to set some specific parameters
643 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
646 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
647 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
648 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
649 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
650 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
651 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
652 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
653 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
654 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
655 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
656 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
657 if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
658 if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
659 if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
660 if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
661 if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
662 if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
663 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
664 if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
665 if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
667 if (config->GetValue("ProcessOld")) fProcessOld = (Bool_t)((TObjString*)config->GetValue("ProcessOld"))->GetString().Atoi();
668 if (config->GetValue("ProcessNew")) fProcessNew = (Bool_t)((TObjString*)config->GetValue("ProcessNew"))->GetString().Atoi();
669 if (config->GetValue("AnalyseNew")) fAnalyseNew = (Bool_t)((TObjString*)config->GetValue("AnalyseNew"))->GetString().Atoi();
671 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
672 for (Int_t i=0;i<5;++i){
678 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
681 //_____________________________________________________________________
682 AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
685 // assignment operator
687 if (&source == this) return *this;
688 new (this) AliTPCCalibCE(source);
692 //_____________________________________________________________________
693 AliTPCCalibCE::~AliTPCCalibCE()
699 fCalRocArrayT0.Delete();
700 fCalRocArrayT0Err.Delete();
701 fCalRocArrayQ.Delete();
702 fCalRocArrayRMS.Delete();
703 fCalRocArrayOutliers.Delete();
705 fHistoQArray.Delete();
706 fHistoT0Array.Delete();
707 fHistoRMSArray.Delete();
709 fHistoTmean.Delete();
711 fParamArrayEventPol1.Delete();
712 fParamArrayEventPol2.Delete();
713 fTMeanArrayEvent.Delete();
714 fQMeanArrayEvent.Delete();
716 fPadTimesArrayEvent.Delete();
717 fPadQArrayEvent.Delete();
718 fPadRMSArrayEvent.Delete();
719 fPadPedestalArrayEvent.Delete();
721 fArrHnDrift.SetOwner();
722 fArrHnDrift.Delete();
725 fArrFitGraphs->SetOwner();
726 delete fArrFitGraphs;
729 //_____________________________________________________________________
730 Int_t AliTPCCalibCE::Update(const Int_t icsector,
733 const Int_t icTimeBin,
734 const Float_t csignal)
737 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
738 // no extra analysis necessary. Assumes knowledge of the signal shape!
739 // assumes that it is looped over consecutive time bins of one pad
742 if (!fProcessOld) return 0;
745 if (icRow<0) return 0;
746 if (icPad<0) return 0;
747 if (icTimeBin<0) return 0;
748 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
750 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
752 //init first pad and sector in this event
753 if ( fCurrentChannel == -1 ) {
755 fCurrentChannel = iChannel;
756 fCurrentSector = icsector;
760 //process last pad if we change to a new one
761 if ( iChannel != fCurrentChannel ){
763 fLastSector=fCurrentSector;
764 fCurrentChannel = iChannel;
765 fCurrentSector = icsector;
769 //fill signals for current pad
770 fPadSignal[icTimeBin]=csignal;
771 if ( csignal > fMaxPadSignal ){
772 fMaxPadSignal = csignal;
773 fMaxTimeBin = icTimeBin;
778 //_____________________________________________________________________
779 void AliTPCCalibCE::ProcessBunch(const Int_t sector, const Int_t row, const Int_t pad,
780 const Int_t length, const UInt_t startTimeBin, const UShort_t* signal)
783 // new filling method to fill the THnSparse histogram
786 //only in new processing mode
787 if (!fProcessNew) return;
788 //don't use the IROCs
789 if (sector<36) return;
790 //only bunches with reasonable length
791 if (length<3||length>10) return;
793 UShort_t timeBin = (UShort_t)startTimeBin;
794 //skip first laser layer
795 if (timeBin<200) return;
796 Double_t timeBurst=SetBurstHnDrift();
798 //after 1 event setup peak ranges
799 if (fNevents==1&&fPeaks[4]==0) {
801 fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
804 fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
807 // After the first event only fill every 5th bin in a row with the CE information
809 if (fPeaks[4]>100&&TMath::Abs((Short_t)fPeaks[4]-(Short_t)timeBin)<(Short_t)fPeakWidths[4]){
816 if (!IsPeakInRange(timeBin+length/2)) return;
818 Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
819 (Double_t)padFill,(Double_t)timeBin,timeBurst};
821 for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
822 Float_t sig=(Float_t)signal[iTimeBin];
824 fHnDrift->Fill(x,sig);
828 //_____________________________________________________________________
829 void AliTPCCalibCE::FindLaserLayers()
832 // Find the laser layer positoins
836 //find CE signal position and width
837 TH1D *hproj=fHnDrift->Projection(3);
838 hproj->GetXaxis()->SetRangeUser(700,1030);
839 Int_t maxbin=hproj->GetMaximumBin();
840 Double_t binc=hproj->GetBinCenter(maxbin);
841 hproj->GetXaxis()->SetRangeUser(binc-10,binc+10);
843 fPeaks[4]=(UShort_t)TMath::Nint(hproj->GetMean());
844 // fPeakWidths[4]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
845 fPeakWidths[4]=(UShort_t)TMath::Nint(10.);
848 Int_t timepos=fPeaks[4]-2*fPeakWidths[4];
851 for (Int_t i=3; i>=0; --i){
852 hproj->GetXaxis()->SetRangeUser(timepos-width,timepos);
853 fPeaks[i]=(UShort_t)TMath::Nint(hproj->GetMean());
854 fPeakWidths[i]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
856 timepos=fPeaks[i]-width/2;
859 //check width and reset peak if >100
860 for (Int_t i=0; i<5; ++i){
861 if (fPeakWidths[i]>100) {
870 //_____________________________________________________________________
871 void AliTPCCalibCE::FindPedestal(Float_t part)
874 // find pedestal and noise for the current pad. Use either database or
875 // truncated mean with part*100%
877 Bool_t noPedestal = kTRUE;
879 //use pedestal database if set
880 if (fPedestalTPC&&fPadNoiseTPC){
881 //only load new pedestals if the sector has changed
882 if ( fCurrentSector!=fLastSector ){
883 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
884 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
887 if ( fPedestalROC&&fPadNoiseROC ){
888 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
889 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
895 //if we are not running with pedestal database, or for the current sector there is no information
896 //available, calculate the pedestal and noise on the fly
900 if ( fIsZeroSuppressed ) return;
901 const Int_t kPedMax = 100; //maximum pedestal value
910 UShort_t histo[kPedMax];
911 memset(histo,0,kPedMax*sizeof(UShort_t));
913 //fill pedestal histogram
914 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
915 padSignal = fPadSignal[i];
916 if (padSignal<=0) continue;
917 if (padSignal>max && i>10) {
921 if (padSignal>kPedMax-1) continue;
922 histo[int(padSignal+0.5)]++;
926 for (Int_t i=1; i<kPedMax; ++i){
927 if (count1<count0*0.5) median=i;
932 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
934 for (Int_t idelta=1; idelta<10; ++idelta){
935 if (median-idelta<=0) continue;
936 if (median+idelta>kPedMax) continue;
937 if (count<part*count1){
938 count+=histo[median-idelta];
939 mean +=histo[median-idelta]*(median-idelta);
940 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
941 count+=histo[median+idelta];
942 mean +=histo[median+idelta]*(median+idelta);
943 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
948 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
954 //_____________________________________________________________________
955 void AliTPCCalibCE::UpdateCETimeRef()
957 // Find the time reference of the last valid CE signal in sector
958 // for irocs of the A-Side the reference of the corresponging OROC is returned
959 // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
960 if ( fLastSector == fCurrentSector ) return;
961 Int_t sector=fCurrentSector;
962 if ( sector < 18 ) sector+=36;
964 TVectorF *vtRef = GetTMeanEvents(sector);
965 if ( !vtRef ) return;
966 Int_t vtRefSize= vtRef->GetNrows();
967 if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
968 else vtRefSize=fNevents;
969 while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
970 fCurrentCETimeRef=(*vtRef)[vtRefSize];
971 AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
973 //_____________________________________________________________________
974 void AliTPCCalibCE::FindCESignal(TVectorD ¶m, Float_t &qSum, const TVectorF maxima)
977 // Find position, signal width and height of the CE signal (last signal)
978 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
979 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
982 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
984 Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise; // threshold for the signal sum
985 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
986 const Int_t kCemax = fPeakIntPlus;
988 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
990 // find maximum closest to the sector mean from the last event
991 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
992 // get sector mean of last event
993 Float_t tmean = fCurrentCETimeRef;
994 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
995 minDist = tmean-maxima[imax];
996 cemaxpos = (Int_t)maxima[imax];
999 // printf("L1 phase TB: %f\n",GetL1PhaseTB());
1001 ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
1002 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
1003 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
1004 Float_t signal = fPadSignal[i]-fPadPedestal;
1006 ceTime+=signal*(i+0.5);
1007 ceRMS +=signal*(i+0.5)*(i+0.5);
1013 if (ceQmax&&ceQsum>ceSumThreshold) {
1015 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
1016 ceTime-=GetL1PhaseTB();
1017 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
1018 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
1020 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
1021 // the pick-up signal should scale with the pad area. In addition
1022 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
1023 // ratio 2/3. The pad area we express in cm2. We normalise the signal
1024 // to the OROC signal (factor 2/3 for the IROCs).
1025 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
1026 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
1029 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
1030 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
1042 //_____________________________________________________________________
1043 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
1046 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
1047 // and 'tplus' timebins after 'pos'
1049 if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
1050 for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
1051 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
1052 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
1053 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
1055 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
1059 //_____________________________________________________________________
1060 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
1063 // Find local maxima on the pad signal and Histogram them
1065 Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal
1068 for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
1069 if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
1070 if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
1071 if (count<maxima.GetNrows()){
1072 maxima.GetMatrixArray()[count++]=i;
1073 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
1074 i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin fPeakDetMinus+fPeakDetPlus-1
1079 //_____________________________________________________________________
1080 void AliTPCCalibCE::ProcessPad()
1083 // Process data of current pad
1087 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
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
1090 FindLocalMaxima(maxima);
1091 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
1093 UpdateCETimeRef(); // update the time refenrence for the current sector
1094 if ( fCurrentCETimeRef<1e-30 ) return; //return if we don't have time 0 info, eg if only one side has laser
1097 FindCESignal(param, qSum, maxima);
1099 Double_t meanT = param[1];
1100 Double_t sigmaT = param[2];
1102 //Fill Event T0 counter
1103 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
1106 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
1108 //Fill RMS histogram
1109 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
1112 //Fill debugging info
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;
1121 //_____________________________________________________________________
1122 void AliTPCCalibCE::EndEvent()
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
1130 if (fProcessNew) ++fNevents;
1134 //check if last pad has allready been processed, if not do so
1135 if ( fMaxTimeBin>-1 ) ProcessPad();
1137 AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
1140 TMatrixD dummy(3,3);
1141 // TVectorF vMeanTime(72);
1142 // TVectorF vMeanQ(72);
1143 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
1144 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
1146 //find mean time0 offset for side A and C
1147 //use only orocs due to the better statistics
1148 Double_t time0Side[2]; //time0 for side A:0 and C:1
1149 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:1
1150 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
1151 for ( Int_t iSec = 36; iSec<72; ++iSec ){
1152 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
1153 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
1155 if ( time0SideCount[0] >0 )
1156 time0Side[0]/=time0SideCount[0];
1157 if ( time0SideCount[1] >0 )
1158 time0Side[1]/=time0SideCount[1];
1159 // end find time0 offset
1160 AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1162 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1163 for ( Int_t iSec = 0; iSec<72; ++iSec ){
1164 AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1165 //find median and then calculate the mean around it
1166 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
1167 if ( !hMeanT ) continue;
1168 //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1169 if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1171 AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1175 Double_t entries = hMeanT->GetEffectiveEntries();
1177 Short_t *arr = hMeanT->GetArray()+1;
1179 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1181 if ( sum>=(entries/2.) ) break;
1184 Int_t firstBin = fFirstTimeBin+ibin-delta;
1185 Int_t lastBin = fFirstTimeBin+ibin+delta;
1186 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1187 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
1188 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1190 // check boundaries for ebye info of mean time
1191 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1192 Int_t vSize=vMeanTime->GetNrows();
1193 if ( vSize < fNevents+1 ){
1194 vMeanTime->ResizeTo(vSize+100);
1197 // store mean time for the readout sides
1198 vSize=fVTime0SideA.GetNrows();
1199 if ( vSize < fNevents+1 ){
1200 fVTime0SideA.ResizeTo(vSize+100);
1201 fVTime0SideC.ResizeTo(vSize+100);
1203 fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1204 fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
1206 vMeanTime->GetMatrixArray()[fNevents]=median;
1210 TVectorF *vTimes = GetPadTimesEvent(iSec);
1211 if ( !vTimes ) continue; //continue if no time information for this sector is available
1213 AliTPCCalROC calIrocOutliers(0);
1214 AliTPCCalROC calOrocOutliers(36);
1216 // calculate mean Q of the sector
1217 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1218 vSize=vMeanQ->GetNrows();
1219 if ( vSize < fNevents+1 ){
1220 vMeanQ->ResizeTo(vSize+100);
1223 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1224 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1226 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1227 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
1229 //set values for temporary roc calibration class
1231 calIroc->SetValue(iChannel, time);
1232 if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
1235 calOroc->SetValue(iChannel, time);
1236 if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
1239 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1240 // test that we really found the CE signal reliably
1241 if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1242 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
1246 //------------------------------- Debug start ------------------------------
1247 if ( GetStreamLevel()>0 ){
1248 TTreeSRedirector *streamer=GetDebugStreamer();
1254 Float_t q = (*GetPadQEvent(iSec))[iChannel];
1255 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1257 UInt_t channel=iChannel;
1260 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1261 pad = channel-fROC->GetRowIndexes(sector)[row];
1262 padc = pad-(fROC->GetNPads(sector,row)/2);
1264 // TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1265 // Form("hSignalD%d.%d.%d",sector,row,pad),
1266 // fLastTimeBin-fFirstTimeBin,
1267 // fFirstTimeBin,fLastTimeBin);
1268 // h1->SetDirectory(0);
1270 // for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1271 // h1->Fill(i,fPadSignal(i));
1274 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1275 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1276 Double_t t0Side = time0Side[(iSec/18)%2];
1277 (*streamer) << "DataPad" <<
1278 "Event=" << fNevents <<
1279 "RunNumber=" << fRunNumber <<
1280 "TimeStamp=" << fTimeStamp <<
1281 "Sector="<< sector <<
1285 "PadSec="<< channel <<
1286 "Time0Sec=" << t0Sec <<
1287 "Time0Side=" << t0Side <<
1291 "MeanQ=" << meanQ <<
1292 // "hist.=" << h1 <<
1298 //----------------------------- Debug end ------------------------------
1299 }// end channel loop
1302 //do fitting now only in debug mode
1303 if (GetDebugLevel()>0){
1304 TVectorD paramPol1(3);
1305 TVectorD paramPol2(6);
1306 TMatrixD matPol1(3,3);
1307 TMatrixD matPol2(6,6);
1311 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1313 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1314 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1316 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1317 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1320 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1321 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1324 //------------------------------- Debug start ------------------------------
1325 if ( GetStreamLevel()>0 ){
1326 TTreeSRedirector *streamer=GetDebugStreamer();
1328 (*streamer) << "DataRoc" <<
1329 // "Event=" << fEvent <<
1330 "RunNumber=" << fRunNumber <<
1331 "TimeStamp=" << fTimeStamp <<
1333 "hMeanT.=" << hMeanT <<
1334 "median=" << median <<
1335 "paramPol1.=" << ¶mPol1 <<
1336 "paramPol2.=" << ¶mPol2 <<
1337 "matPol1.=" << &matPol1 <<
1338 "matPol2.=" << &matPol2 <<
1339 "chi2Pol1=" << chi2Pol1 <<
1340 "chi2Pol2=" << chi2Pol2 <<
1345 //------------------------------- Debug end ------------------------------
1348 //return if no sector has a valid mean time
1349 if ( nSecMeanT == 0 ) return;
1352 // fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1353 // fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1354 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1355 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1356 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1358 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1359 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1362 fOldRunNumber = fRunNumber;
1366 AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1368 //_____________________________________________________________________
1369 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1370 Int_t nbinsY, Float_t ymin, Float_t ymax,
1371 const Char_t *type, Bool_t force)
1374 // return pointer to TH2S histogram of 'type'
1375 // if force is true create a new histogram if it doesn't exist allready
1377 if ( !force || arr->UncheckedAt(sector) )
1378 return (TH2S*)arr->UncheckedAt(sector);
1380 // if we are forced and histogram doesn't exist yet create it
1381 // new histogram with Q calib information. One value for each pad!
1382 TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1384 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1385 hist->SetDirectory(0);
1386 arr->AddAt(hist,sector);
1389 //_____________________________________________________________________
1390 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1393 // return pointer to T0 histogram
1394 // if force is true create a new histogram if it doesn't exist allready
1396 TObjArray *arr = &fHistoT0Array;
1397 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1399 //_____________________________________________________________________
1400 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1403 // return pointer to Q histogram
1404 // if force is true create a new histogram if it doesn't exist allready
1406 TObjArray *arr = &fHistoQArray;
1407 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1409 //_____________________________________________________________________
1410 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1413 // return pointer to Q histogram
1414 // if force is true create a new histogram if it doesn't exist allready
1416 TObjArray *arr = &fHistoRMSArray;
1417 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1419 //_____________________________________________________________________
1420 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1421 const Char_t *type, Bool_t force)
1424 // return pointer to TH1S histogram
1425 // if force is true create a new histogram if it doesn't exist allready
1427 if ( !force || arr->UncheckedAt(sector) )
1428 return (TH1S*)arr->UncheckedAt(sector);
1430 // if we are forced and histogram doesn't yes exist create it
1431 // new histogram with calib information. One value for each pad!
1432 TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1433 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1434 hist->SetDirectory(0);
1435 arr->AddAt(hist,sector);
1438 //_____________________________________________________________________
1439 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1442 // return pointer to Q histogram
1443 // if force is true create a new histogram if it doesn't exist allready
1445 TObjArray *arr = &fHistoTmean;
1446 return GetHisto(sector, arr, "LastTmean", force);
1448 //_____________________________________________________________________
1449 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1452 // return pointer to Pad Info from 'arr' for the current event and sector
1453 // if force is true create it if it doesn't exist allready
1455 if ( !force || arr->UncheckedAt(sector) )
1456 return (TVectorF*)arr->UncheckedAt(sector);
1458 TVectorF *vect = new TVectorF(size);
1459 arr->AddAt(vect,sector);
1462 //_____________________________________________________________________
1463 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1466 // return pointer to Pad Times Array for the current event and sector
1467 // if force is true create it if it doesn't exist allready
1469 TObjArray *arr = &fPadTimesArrayEvent;
1470 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1472 //_____________________________________________________________________
1473 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1476 // return pointer to Pad Q Array for the current event and sector
1477 // if force is true create it if it doesn't exist allready
1478 // for debugging purposes only
1481 TObjArray *arr = &fPadQArrayEvent;
1482 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1484 //_____________________________________________________________________
1485 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1488 // return pointer to Pad RMS Array for the current event and sector
1489 // if force is true create it if it doesn't exist allready
1490 // for debugging purposes only
1492 TObjArray *arr = &fPadRMSArrayEvent;
1493 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1495 //_____________________________________________________________________
1496 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1499 // return pointer to Pad RMS Array for the current event and sector
1500 // if force is true create it if it doesn't exist allready
1501 // for debugging purposes only
1503 TObjArray *arr = &fPadPedestalArrayEvent;
1504 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1506 //_____________________________________________________________________
1507 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1510 // return pointer to the EbyE info of the mean arrival time for 'sector'
1511 // if force is true create it if it doesn't exist allready
1513 TObjArray *arr = &fTMeanArrayEvent;
1514 return GetVectSector(sector,arr,100,force);
1516 //_____________________________________________________________________
1517 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1520 // return pointer to the EbyE info of the mean arrival time for 'sector'
1521 // if force is true create it if it doesn't exist allready
1523 TObjArray *arr = &fQMeanArrayEvent;
1524 return GetVectSector(sector,arr,100,force);
1526 //_____________________________________________________________________
1527 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1530 // return pointer to ROC Calibration
1531 // if force is true create a new histogram if it doesn't exist allready
1533 if ( !force || arr->UncheckedAt(sector) )
1534 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1536 // if we are forced and histogram doesn't yes exist create it
1538 // new AliTPCCalROC for T0 information. One value for each pad!
1539 AliTPCCalROC *croc = new AliTPCCalROC(sector);
1540 arr->AddAt(croc,sector);
1543 //_____________________________________________________________________
1544 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1547 // return pointer to Time 0 ROC Calibration
1548 // if force is true create a new histogram if it doesn't exist allready
1550 TObjArray *arr = &fCalRocArrayT0;
1551 return GetCalRoc(sector, arr, force);
1553 //_____________________________________________________________________
1554 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1557 // return pointer to the error of Time 0 ROC Calibration
1558 // if force is true create a new histogram if it doesn't exist allready
1560 TObjArray *arr = &fCalRocArrayT0Err;
1561 return GetCalRoc(sector, arr, force);
1563 //_____________________________________________________________________
1564 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1567 // return pointer to T0 ROC Calibration
1568 // if force is true create a new histogram if it doesn't exist allready
1570 TObjArray *arr = &fCalRocArrayQ;
1571 return GetCalRoc(sector, arr, force);
1573 //_____________________________________________________________________
1574 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1577 // return pointer to signal width ROC Calibration
1578 // if force is true create a new histogram if it doesn't exist allready
1580 TObjArray *arr = &fCalRocArrayRMS;
1581 return GetCalRoc(sector, arr, force);
1583 //_____________________________________________________________________
1584 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1587 // return pointer to Outliers
1588 // if force is true create a new histogram if it doesn't exist allready
1590 TObjArray *arr = &fCalRocArrayOutliers;
1591 return GetCalRoc(sector, arr, force);
1593 //_____________________________________________________________________
1594 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1597 // return pointer to TObjArray of fit parameters
1598 // if force is true create a new histogram if it doesn't exist allready
1600 if ( !force || arr->UncheckedAt(sector) )
1601 return (TObjArray*)arr->UncheckedAt(sector);
1603 // if we are forced and array doesn't yes exist create it
1605 // new TObjArray for parameters
1606 TObjArray *newArr = new TObjArray;
1607 arr->AddAt(newArr,sector);
1610 //_____________________________________________________________________
1611 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1614 // return pointer to TObjArray of fit parameters from plane fit
1615 // if force is true create a new histogram if it doesn't exist allready
1617 TObjArray *arr = &fParamArrayEventPol1;
1618 return GetParamArray(sector, arr, force);
1620 //_____________________________________________________________________
1621 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1624 // return pointer to TObjArray of fit parameters from parabola fit
1625 // if force is true create a new histogram if it doesn't exist allready
1627 TObjArray *arr = &fParamArrayEventPol2;
1628 return GetParamArray(sector, arr, force);
1631 //_____________________________________________________________________
1632 void AliTPCCalibCE::CreateDVhist()
1635 // Setup the THnSparse for the drift velocity determination
1639 //roc, row, pad, timebin, timestamp
1640 TTimeStamp begin(2010,01,01,0,0,0);
1641 TTimeStamp end(2035,01,01,0,0,0);
1642 Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1644 Int_t bins[kHnBinsDV] = { 72, 96, 140, 1030, nbinsTime};
1645 Double_t xmin[kHnBinsDV] = { 0., 0., 0., 0., (Double_t)begin.GetSec()};
1646 Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1648 fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1649 fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1650 fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1651 fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1652 fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1653 fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
1657 //_____________________________________________________________________
1658 void AliTPCCalibCE::ResetEvent()
1661 // Reset global counters -- Should be called before each event is processed
1670 fPadTimesArrayEvent.Delete();
1671 fPadQArrayEvent.Delete();
1672 fPadRMSArrayEvent.Delete();
1673 fPadPedestalArrayEvent.Delete();
1675 for ( Int_t i=0; i<72; ++i ){
1676 fVTime0Offset.GetMatrixArray()[i]=0;
1677 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1678 fVMeanQ.GetMatrixArray()[i]=0;
1679 fVMeanQCounter.GetMatrixArray()[i]=0;
1682 //_____________________________________________________________________
1683 void AliTPCCalibCE::ResetPad()
1686 // Reset pad infos -- Should be called after a pad has been processed
1688 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1695 //_____________________________________________________________________
1696 void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
1699 // Merge ce to the current AliTPCCalibCE
1702 Int_t nCEevents = ce->GetNeventsProcessed();
1704 if (fProcessOld&&ce->fProcessOld){
1706 for (Int_t iSec=0; iSec<72; ++iSec){
1707 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1708 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1709 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1713 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1714 TH2S *hRefQ = GetHistoQ(iSec);
1715 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1717 TH2S *hist = new TH2S(*hRefQmerge);
1718 hist->SetDirectory(0);
1719 fHistoQArray.AddAt(hist, iSec);
1721 hRefQmerge->SetDirectory(dir);
1724 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1725 TH2S *hRefT0 = GetHistoT0(iSec);
1726 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1728 TH2S *hist = new TH2S(*hRefT0merge);
1729 hist->SetDirectory(0);
1730 fHistoT0Array.AddAt(hist, iSec);
1732 hRefT0merge->SetDirectory(dir);
1734 if ( hRefRMSmerge ){
1735 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1736 TH2S *hRefRMS = GetHistoRMS(iSec);
1737 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1739 TH2S *hist = new TH2S(*hRefRMSmerge);
1740 hist->SetDirectory(0);
1741 fHistoRMSArray.AddAt(hist, iSec);
1743 hRefRMSmerge->SetDirectory(dir);
1748 // merge time information
1751 for (Int_t iSec=0; iSec<72; ++iSec){
1752 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1753 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1754 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1755 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1757 TObjArray *arrPol1 = 0x0;
1758 TObjArray *arrPol2 = 0x0;
1759 TVectorF *vMeanTime = 0x0;
1760 TVectorF *vMeanQ = 0x0;
1763 if ( arrPol1CE && arrPol2CE ){
1764 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1765 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1766 arrPol1->Expand(fNevents+nCEevents);
1767 arrPol2->Expand(fNevents+nCEevents);
1769 if ( vMeanTimeCE && vMeanQCE ){
1770 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1771 vMeanQ = GetQMeanEvents(iSec,kTRUE);
1772 vMeanTime->ResizeTo(fNevents+nCEevents);
1773 vMeanQ->ResizeTo(fNevents+nCEevents);
1776 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1777 if ( arrPol1CE && arrPol2CE ){
1778 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1779 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1780 if ( paramPol1 && paramPol2 ){
1781 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1782 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1785 if ( vMeanTimeCE && vMeanQCE ){
1786 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1787 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1794 const TVectorD& eventTimes = ce->fVEventTime;
1795 const TVectorD& eventIds = ce->fVEventNumber;
1796 const TVectorF& time0SideA = ce->fVTime0SideA;
1797 const TVectorF& time0SideC = ce->fVTime0SideC;
1798 fVEventTime.ResizeTo(fNevents+nCEevents);
1799 fVEventNumber.ResizeTo(fNevents+nCEevents);
1800 fVTime0SideA.ResizeTo(fNevents+nCEevents);
1801 fVTime0SideC.ResizeTo(fNevents+nCEevents);
1803 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1804 Double_t evTime = eventTimes.GetMatrixArray()[iEvent];
1805 Double_t evId = eventIds.GetMatrixArray()[iEvent];
1806 Float_t t0SideA = time0SideA.GetMatrixArray()[iEvent];
1807 Float_t t0SideC = time0SideC.GetMatrixArray()[iEvent];
1809 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1810 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1811 fVTime0SideA.GetMatrixArray()[fNevents+iEvent] = t0SideA;
1812 fVTime0SideC.GetMatrixArray()[fNevents+iEvent] = t0SideC;
1816 if (fProcessNew&&ce->fProcessNew) {
1817 if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1818 AliError("Number of bursts in the instances to merge are different. No merging done!");
1820 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1821 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1822 THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1823 if (h && hce) h->Add(hce);
1824 else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
1826 //TODO: What to do with fTimeBursts???
1830 fNevents+=nCEevents; //increase event counter
1833 //_____________________________________________________________________
1834 Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1837 // Merge all objects of this type in list
1843 AliTPCCalibCE *ce=0;
1846 while ( (o=next()) ){
1847 ce=dynamic_cast<AliTPCCalibCE*>(o);
1857 //_____________________________________________________________________
1858 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1861 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1862 // or side (-1: A-Side, -2: C-Side)
1863 // xVariable: 0-event time, 1-event id, 2-internal event counter
1864 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1865 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1866 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1867 // not used for mean time and mean Q )
1868 // for an example see class description at the beginning
1871 TVectorD *xVar = 0x0;
1872 TObjArray *aType = 0x0;
1876 if ( (sector<-2) || (sector>71) ) return 0x0; //sector outside valid range
1877 if ( (xVariable<0) || (xVariable>2) ) return 0x0; //invalid x-variable
1878 if ( (fitType<0) || (fitType>3) ) return 0x0; //invalid fit type
1879 if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1880 if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1881 if ( sector<0 && fitType!=2) return 0x0; //for side wise information only mean time is available
1885 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1886 aType = &fParamArrayEventPol1;
1887 if ( aType->At(sector)==0x0 ) return 0x0;
1889 else if ( fitType==1 ){
1890 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1891 aType = &fParamArrayEventPol2;
1892 if ( aType->At(sector)==0x0 ) return 0x0;
1896 if ( xVariable == 0 ) xVar = &fVEventTime;
1897 if ( xVariable == 1 ) xVar = &fVEventNumber;
1898 if ( xVariable == 2 ) {
1899 xVar = new TVectorD(fNevents);
1900 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1903 Double_t *x = new Double_t[fNevents];
1904 Double_t *y = new Double_t[fNevents];
1906 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1908 TObjArray *events = (TObjArray*)(aType->At(sector));
1909 if ( events->GetSize()<=ievent ) break;
1910 TVectorD *v = (TVectorD*)(events->At(ievent));
1911 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1912 } else if (fitType == 2) {
1913 Double_t xValue=(*xVar)[ievent];
1915 if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1916 else if (sector==-1) yValue=fVTime0SideA(ievent);
1917 else if (sector==-2) yValue=fVTime0SideC(ievent);
1918 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1919 }else if (fitType == 3) {
1920 Double_t xValue=(*xVar)[ievent];
1921 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1922 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1926 TGraph *gr = new TGraph(npoints);
1927 //sort xVariable increasing
1928 Int_t *sortIndex = new Int_t[npoints];
1929 TMath::Sort(npoints,x,sortIndex, kFALSE);
1930 for (Int_t i=0;i<npoints;++i){
1931 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1935 if ( xVariable == 2 ) delete xVar;
1938 delete [] sortIndex;
1941 //_____________________________________________________________________
1942 void AliTPCCalibCE::Analyse()
1945 // Calculate calibration constants
1950 TVectorD paramT0(3);
1951 TVectorD paramRMS(3);
1952 TMatrixD dummy(3,3);
1954 Float_t channelCounter=0;
1959 for (Int_t iSec=0; iSec<72; ++iSec){
1960 TH2S *hT0 = GetHistoT0(iSec);
1961 if (!hT0 ) continue;
1963 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1964 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1965 AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1966 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1967 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1969 TH2S *hQ = GetHistoQ(iSec);
1970 TH2S *hRMS = GetHistoRMS(iSec);
1972 Short_t *arrayhQ = hQ->GetArray();
1973 Short_t *arrayhT0 = hT0->GetArray();
1974 Short_t *arrayhRMS = hRMS->GetArray();
1976 UInt_t nChannels = fROC->GetNChannels(iSec);
1984 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1987 Float_t cogTime0 = -1000;
1988 Float_t cogQ = -1000;
1989 Float_t cogRMS = -1000;
1995 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1996 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1997 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1999 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2001 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2003 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2008 //outlier specifications
2009 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2016 rocQ->SetValue(iChannel, cogQ*cogQ);
2017 rocT0->SetValue(iChannel, cogTime0);
2018 rocT0Err->SetValue(iChannel, rmsT0);
2019 rocRMS->SetValue(iChannel, cogRMS);
2020 rocOut->SetValue(iChannel, cogOut);
2024 if ( GetStreamLevel() > 0 ){
2025 TTreeSRedirector *streamer=GetDebugStreamer();
2028 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2029 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2030 padc = pad-(fROC->GetNPads(iSec,row)/2);
2032 (*streamer) << "DataEnd" <<
2033 "Sector=" << iSec <<
2037 "PadSec=" << iChannel <<
2039 "T0=" << cogTime0 <<
2049 if ( channelCounter>0 ){
2050 fMeanT0rms/=channelCounter;
2051 fMeanQrms/=channelCounter;
2052 fMeanRMSrms/=channelCounter;
2054 // if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2055 // delete fDebugStreamer;
2056 // fDebugStreamer = 0x0;
2057 fVEventTime.ResizeTo(fNevents);
2058 fVEventNumber.ResizeTo(fNevents);
2059 fVTime0SideA.ResizeTo(fNevents);
2060 fVTime0SideC.ResizeTo(fNevents);
2063 if (fProcessNew&&fAnalyseNew){
2065 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2066 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2067 h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2076 // New functions that also use the laser tracks
2081 //_____________________________________________________________________
2082 void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2085 //Find the local maximums for the projections to each axis
2088 //find laser layer positoins
2089 fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2091 THnSparse *hProj=fHnDrift;
2092 Double_t posCE[4]={0.,0.,0.,0.};
2093 Double_t widthCE[4]={0.,0.,0.,0.};
2095 // if(fPeaks[4]!=0){
2096 // find central electrode position once more, separately for IROC, OROC, A-, C-Side
2098 for (Int_t i=0; i<4; ++i){
2099 hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2100 TH1 *h=fHnDrift->Projection(3);
2101 h->GetXaxis()->SetRangeUser(fPeaks[4]-fPeakWidths[4],fPeaks[4]+fPeakWidths[4]);
2102 Int_t nbinMax=h->GetMaximumBin();
2103 Double_t maximum=h->GetMaximum();
2104 // Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2105 // if (nbinMax<700||maximum<maxExpected) continue;
2106 Double_t xbinMax=h->GetBinCenter(nbinMax);
2107 TF1 fgaus("gaus","gaus",xbinMax-10,xbinMax+10);
2108 fgaus.SetParameters(maximum,xbinMax,2);
2109 fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2110 fgaus.SetParLimits(2,0.2,4.);
2111 h->Fit(&fgaus,"RQN");
2112 // Double_t deltaX=4*fgaus.GetParameter(2);
2113 // xbinMax=fgaus.GetParameter(1);
2115 posCE[i]=fgaus.GetParameter(1);
2116 widthCE[i]=4*fgaus.GetParameter(2);
2117 hProj->GetAxis(0)->SetRangeUser(0,72);
2120 //Current drift velocity
2121 Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2122 // cout<<"vdrift="<<vdrift<<endl;
2124 AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2125 //loop over all entries in the histogram
2127 for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2128 //get entry position and content
2129 Double_t adc=hProj->GetBinContent(ichunk,coord);
2132 Int_t sector = coord[0]-1;
2133 Int_t row = coord[1]-1;
2134 Int_t pad = coord[2]-1;
2135 Int_t timeBin= coord[3]-1;
2136 Double_t time = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2137 Int_t side = (sector/18)%2;
2139 // fPeaks[4]=(UInt_t)posCE[sector/18];
2140 // fPeakWidths[4]=(UInt_t)widthCE[sector/18];
2143 if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2144 if (adc < 5 ) continue;
2145 if (IsEdgePad(sector,row,pad)) continue;
2146 // if (!IsPeakInRange(timeBin)) continue;
2147 // if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2148 // if (isCE&&(sector!=0)) continue;
2150 Int_t padmin=-2, padmax=2;
2151 Int_t timemin=-2, timemax=2;
2152 Int_t minsumperpad=2;
2153 //CE or laser tracks
2155 if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2164 // Find local maximum and cogs
2166 Bool_t isMaximum=kTRUE;
2167 Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2168 Double_t cogY=0, rmsY=0;
2171 // for position calculation use
2172 for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2174 fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2176 for(Int_t itime=timemin;itime<=timemax;++itime){
2178 Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2179 Double_t val=hProj->GetBinContent(a);
2182 tbcm +=(timeBin+itime)*val;
2183 padcm+=(pad+ipad)*val;
2186 rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2187 rmspad+=(pad+ipad)*(pad+ipad)*val;
2188 rmsY +=lxyz[1]*lxyz[1]*val;
2196 if (!isMaximum) break;
2199 if (!isMaximum||!npart) continue;
2200 if (totalmass<npart*minsumperpad) continue;
2201 if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2202 //for CE we want only one pad by construction
2212 rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2213 rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2214 rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2216 Int_t cog=TMath::Nint(padcm);
2218 // timebin --> z position
2219 Float_t zlength=fParam->GetZLength(side);
2220 // Float_t timePos=tbcm+(1000-fPeaks[4])
2221 // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2222 Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2224 // local to global transformation--> x and y positions
2226 fROC->GetPositionLocal(sector,row,pad,padlxyz);
2228 Double_t gxyz[3]={padlxyz[0],cogY,gz};
2229 Double_t lxyz[3]={padlxyz[0],cogY,gz};
2230 Double_t igxyz[3]={0,0,0};
2232 t1.RotatedGlobal2Global(sector,gxyz);
2238 //find track id and index of the position in the track (row)
2241 index=row+(sector>35)*fROC->GetNRows(0);
2242 trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2244 trackID=336+((sector/18)%2);
2245 index= fROC->GetRowIndexes(sector)[row]+pad; // global pad position in sector
2247 index+=(sector%18)*fROC->GetNChannels(sector);
2249 index+=18*fROC->GetNChannels(0);
2250 index+=(sector%18)*fROC->GetNChannels(sector);
2252 //TODO: find out about the multiple peaks in the CE
2253 // mindist=TMath::Abs(fPeaks[4]-tbcm);
2257 // fill track vectors
2259 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2260 Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2263 // travel time effect of light includes
2265 Double_t raylength=ltr->GetRayLength();
2266 Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2267 ltr->GetXYZ(globmir);
2270 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2271 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2272 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2275 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2276 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2277 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2281 if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2282 ltr->fVecSec->GetMatrixArray()[index]=sector;
2283 ltr->fVecP2->GetMatrixArray()[index]=0;
2284 ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2285 ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2286 ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2287 ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2288 ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2289 ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2290 ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2291 // ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2293 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2294 ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2295 igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2296 igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2297 igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2301 if (fStreamLevel>4){
2302 (*GetDebugStreamer()) << "clusters" <<
2303 "run=" << fRunNumber <<
2304 "timestamp=" << timestamp <<
2305 "burst=" << burst <<
2311 "timebin=" << timeBin <<
2312 "cogCE=" << posCE[sector/18] <<
2313 "withCE=" << widthCE[sector/18] <<
2314 "index=" << index <<
2316 "padcm=" << padcm <<
2317 "rmspad=" << rmspad <<
2320 "rmstb=" << rmstb <<
2324 "lx=" << padlxyz[0]<<
2326 "lypad=" << padlxyz[1]<<
2333 "igx=" << igxyz[0] <<
2334 "igy=" << igxyz[1] <<
2335 "igz=" << igxyz[2] <<
2337 "mind=" << mindist <<
2339 "trackid=" << trackID <<
2340 "trackid2=" << trackID2 <<
2341 "npart=" << npart <<
2343 } // end stream levelmgz.fElements
2349 //_____________________________________________________________________
2350 void AliTPCCalibCE::AnalyseTrack()
2353 // Analyse the tracks
2357 AliTPCLaserTrack::LoadTracks();
2358 // AliTPCParam *param=0x0;
2360 // AliCDBManager *man=AliCDBManager::Instance();
2361 // if (man->GetDefaultStorage()){
2362 // AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2364 // entry->SetOwner(kTRUE);
2365 // param = (AliTPCParam*)(entry->GetObject()->Clone());
2369 // if (fParam) delete fParam;
2372 // AliError("Could not get updated AliTPCParam from OCDB!!!");
2375 //Measured and ideal laser tracks
2376 TObjArray* arrMeasured = SetupMeasured();
2377 TObjArray* arrIdeal = AliTPCLaserTrack::GetTracks();
2378 AddCEtoIdeal(arrIdeal);
2380 //find bursts and loop over them
2381 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2382 Double_t timestamp=fTimeBursts[iburst];
2383 AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2384 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2385 if (!fHnDrift) continue;
2386 UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2387 if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2388 fBinsLastAna[iburst]=entries;
2390 for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2391 // if (iburst==0) FindLaserLayers();
2393 //reset laser tracks
2394 ResetMeasured(arrMeasured);
2396 //find clusters and associate to the tracks
2397 FindLocalMaxima(arrMeasured, timestamp, iburst);
2399 //calculate drift velocity
2400 CalculateDV(arrIdeal,arrMeasured,iburst);
2402 //Dump information to file if requested
2403 if (fStreamLevel>2){
2404 printf("make tree\n");
2405 //laser track information
2407 for (Int_t itrack=0; itrack<338; ++itrack){
2408 TObject *iltr=arrIdeal->UncheckedAt(itrack);
2409 TObject *mltr=arrMeasured->UncheckedAt(itrack);
2410 (*GetDebugStreamer()) << "tracks" <<
2411 "run=" << fRunNumber <<
2412 "time=" << timestamp <<
2413 "burst="<< iburst <<
2420 if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2423 //_____________________________________________________________________
2424 Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2425 const Double_t *peakposloc, Int_t &itrackMin2)
2428 // Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2432 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2433 TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2438 Int_t lastbeam=336/2;
2439 if ( (sector/18)%2 ) {
2446 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2447 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2448 // if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2449 vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2450 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2451 if (TMath::Abs(deltaZ)>40) continue;
2452 vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2453 vSt.RotateZ(ltr->GetAlpha());
2454 vDir.RotateZ(ltr->GetAlpha());
2456 Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2465 Float_t mindist2=10;
2466 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2467 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2468 if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2470 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2471 if (TMath::Abs(deltaZ)>40) continue;
2473 Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2474 if (dist>1) continue;
2486 //_____________________________________________________________________
2487 Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2490 // return true if pad is on the edge of a row
2493 if ( pad == edge1 ) return kTRUE;
2494 Int_t edge2 = fROC->GetNPads(sector,row)-1;
2495 if ( pad == edge2 ) return kTRUE;
2500 //_____________________________________________________________________
2501 TObjArray* AliTPCCalibCE::SetupMeasured()
2504 // setup array of measured laser tracks and CE
2507 TObjArray *arrIdeal = AliTPCLaserTrack::GetTracks();
2508 TObjArray *arrMeasured = new TObjArray(338);
2509 arrMeasured->SetOwner();
2510 for(Int_t itrack=0;itrack<336;itrack++){
2511 arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2515 AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2518 ltrce->fVecSec=new TVectorD(557568/2);
2519 ltrce->fVecP2=new TVectorD(557568/2);
2520 ltrce->fVecPhi=new TVectorD(557568/2);
2521 ltrce->fVecGX=new TVectorD(557568/2);
2522 ltrce->fVecGY=new TVectorD(557568/2);
2523 ltrce->fVecGZ=new TVectorD(557568/2);
2524 ltrce->fVecLX=new TVectorD(557568/2);
2525 ltrce->fVecLY=new TVectorD(557568/2);
2526 ltrce->fVecLZ=new TVectorD(557568/2);
2528 arrMeasured->AddAt(ltrce,336); //CE A-Side
2530 ltrce=new AliTPCLaserTrack;
2533 ltrce->fVecSec=new TVectorD(557568/2);
2534 ltrce->fVecP2=new TVectorD(557568/2);
2535 ltrce->fVecPhi=new TVectorD(557568/2);
2536 ltrce->fVecGX=new TVectorD(557568/2);
2537 ltrce->fVecGY=new TVectorD(557568/2);
2538 ltrce->fVecGZ=new TVectorD(557568/2);
2539 ltrce->fVecLX=new TVectorD(557568/2);
2540 ltrce->fVecLY=new TVectorD(557568/2);
2541 ltrce->fVecLZ=new TVectorD(557568/2);
2542 arrMeasured->AddAt(ltrce,337); //CE C-Side
2547 //_____________________________________________________________________
2548 void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2551 // reset array of measured laser tracks and CE
2553 for(Int_t itrack=0;itrack<338;itrack++){
2554 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2555 ltr->fVecSec->Zero();
2556 ltr->fVecP2->Zero();
2557 ltr->fVecPhi->Zero();
2558 ltr->fVecGX->Zero();
2559 ltr->fVecGY->Zero();
2560 ltr->fVecGZ->Zero();
2561 ltr->fVecLX->Zero();
2562 ltr->fVecLY->Zero();
2563 ltr->fVecLZ->Zero();
2567 //_____________________________________________________________________
2568 void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2571 // Add ideal CE positions to the ideal track data
2576 AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2579 ltrceA->fVecSec=new TVectorD(557568/2);
2580 ltrceA->fVecP2=new TVectorD(557568/2);
2581 ltrceA->fVecPhi=new TVectorD(557568/2);
2582 ltrceA->fVecGX=new TVectorD(557568/2);
2583 ltrceA->fVecGY=new TVectorD(557568/2);
2584 ltrceA->fVecGZ=new TVectorD(557568/2);
2585 ltrceA->fVecLX=new TVectorD(557568/2);
2586 ltrceA->fVecLY=new TVectorD(557568/2);
2587 ltrceA->fVecLZ=new TVectorD(557568/2);
2588 arr->AddAt(ltrceA,336); //CE A-Side
2590 AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2593 ltrceC->fVecSec=new TVectorD(557568/2);
2594 ltrceC->fVecP2=new TVectorD(557568/2);
2595 ltrceC->fVecPhi=new TVectorD(557568/2);
2596 ltrceC->fVecGX=new TVectorD(557568/2);
2597 ltrceC->fVecGY=new TVectorD(557568/2);
2598 ltrceC->fVecGZ=new TVectorD(557568/2);
2599 ltrceC->fVecLX=new TVectorD(557568/2);
2600 ltrceC->fVecLY=new TVectorD(557568/2);
2601 ltrceC->fVecLZ=new TVectorD(557568/2);
2602 arr->AddAt(ltrceC,337); //CE C-Side
2604 //Calculate ideal positoins
2607 Int_t channelSideA=0;
2608 Int_t channelSideC=0;
2609 Int_t channelSide=0;
2610 AliTPCLaserTrack *ltrce=0x0;
2612 for (Int_t isector=0; isector<72; ++isector){
2613 Int_t side=((isector/18)%2);
2614 for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2615 for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2616 fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2617 fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2620 channelSide=channelSideA;
2623 channelSide=channelSideC;
2626 ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2627 ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2628 ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2629 ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2630 ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2631 // ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2632 ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2633 ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2634 // ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2637 ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2638 ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2642 ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2643 ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2653 //_____________________________________________________________________
2654 void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2657 // calculate the drift velocity from the reconstructed clusters associated
2658 // to the ideal laser tracks
2659 // use two different fit scenarios: Separate fit for A- and C-Side
2660 // Common fit for A- and C-Side
2663 if (!fArrFitGraphs){
2664 fArrFitGraphs=new TObjArray;
2667 // static TLinearFitter fdriftA(5,"hyp4");
2668 // static TLinearFitter fdriftC(5,"hyp4");
2669 // static TLinearFitter fdriftAC(6,"hyp5");
2670 Double_t timestamp=fTimeBursts[burst];
2672 static TLinearFitter fdriftA(4,"hyp3");
2673 static TLinearFitter fdriftC(4,"hyp3");
2674 static TLinearFitter fdriftAC(5,"hyp4");
2675 TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2679 Float_t chi2AC = 10;
2684 Double_t minres[3]={20.,1,0.8};
2686 for(Int_t i=0;i<3;i++){
2688 fdriftA.ClearPoints();
2689 fdriftC.ClearPoints();
2690 fdriftAC.ClearPoints();
2699 for (Int_t itrack=0; itrack<338; ++itrack){
2700 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2701 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2703 //-- Exclude the tracks which has the biggest inclanation angle
2704 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2705 Int_t clustercounter=0;
2708 //-- exclude the low intensity tracks
2710 for (Int_t index=0; index<indexMax; ++index){
2712 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2713 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2714 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2716 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2718 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2723 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2725 if (itrack>335) indexMax=557568/2;
2726 for (Int_t index=0; index<indexMax; ++index){
2727 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2728 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2729 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2730 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2732 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2733 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2734 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2735 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2737 //cut if no track info available
2738 if (iltr->GetBundle()==0) continue;
2739 if (iR<133||mR<133) continue;
2740 if(mltr->fVecP2->GetMatrixArray()[index]>minres[i]) continue;
2742 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2743 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2745 //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2746 Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2748 if (iltr->GetSide()==0){
2749 fdriftA.AddPoint(xxx,mdrift,1);
2751 fdriftC.AddPoint(xxx,mdrift,1);
2753 // Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2754 Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2755 fdriftAC.AddPoint(xxx2,mdrift,1);
2758 }//end laser track loop
2768 fdriftA.GetParameters(fitA);
2769 fdriftC.GetParameters(fitC);
2770 fdriftAC.GetParameters(fitAC);
2772 //Parameters: 0 linear offset
2773 // 1 mean drift velocity correction factor
2774 // 2 relative global y gradient
2775 // 3 radial deformation
2776 // 4 IROC/OROC offset
2778 // FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2780 for (Int_t itrack=0; itrack<338; ++itrack){
2781 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2782 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2784 //-- Exclude the tracks which has the biggest inclanation angle
2785 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2786 Int_t clustercounter=0;
2789 //-- exclude the low intensity tracks
2791 for (Int_t index=0; index<indexMax; ++index){
2792 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2793 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2794 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2795 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2797 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2801 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2803 if (itrack>335) indexMax=557568/2;
2804 for (Int_t index=0; index<indexMax; ++index){
2805 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2806 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2807 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2808 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2810 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2811 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2812 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2813 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2815 //cut if no track info available
2816 if (iR<60||mR<60) continue;
2818 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2819 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2822 if (iltr->GetSide()==1) v=&fitC;
2823 // 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);
2824 Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2826 mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2835 //set statistics values
2837 npointsA= fdriftA.GetNpoints();
2838 if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2842 npointsC= fdriftC.GetNpoints();
2843 if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2847 npointsAC= fdriftAC.GetNpoints();
2848 if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2852 if (fStreamLevel>2){
2853 //laser track information
2854 (*GetDebugStreamer()) << "DriftV" <<
2856 "run=" << fRunNumber <<
2857 "time=" << timestamp <<
2858 "fitA.=" << &fitA <<
2859 "fitC.=" << &fitC <<
2860 "fitAC.=" << &fitAC <<
2870 //Parameters: 0 linear offset (global)
2871 // 1 mean drift velocity correction factor
2872 // 2 relative global y gradient
2873 // 3 radial deformation
2874 // 4 IROC/OROC offset
2875 // 5 linear offset relative A-C
2878 TGraphErrors *grA[7];
2879 TGraphErrors *grC[7];
2880 TGraphErrors *grAC[8];
2881 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_");
2882 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_");
2884 TObjArray *arrNames=names.Tokenize(";");
2885 TObjArray *arrNamesAC=namesAC.Tokenize(";");
2888 for (Int_t i=0; i<7; ++i){
2889 TString grName=arrNames->UncheckedAt(i)->GetName();
2891 grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2893 grA[i]=new TGraphErrors;
2894 grA[i]->SetName(grName.Data());
2895 grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2896 fArrFitGraphs->Add(grA[i]);
2898 // Int_t ipoint=grA[i]->GetN();
2900 grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2901 grA[i]->SetPointError(ipoint,60,.0001);
2902 if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2906 for (Int_t i=0; i<7; ++i){
2907 TString grName=arrNames->UncheckedAt(i)->GetName();
2909 grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2911 grC[i]=new TGraphErrors;
2912 grC[i]->SetName(grName.Data());
2913 grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2914 fArrFitGraphs->Add(grC[i]);
2916 // Int_t ipoint=grC[i]->GetN();
2918 grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2919 grC[i]->SetPointError(ipoint,60,.0001);
2920 if (i<4) grC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2924 for (Int_t i=0; i<8; ++i){
2925 TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2927 grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2929 grAC[i]=new TGraphErrors;
2930 grAC[i]->SetName(grName.Data());
2931 grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2932 fArrFitGraphs->Add(grAC[i]);
2934 // Int_t ipoint=grAC[i]->GetN();
2936 grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2937 grAC[i]->SetPointError(ipoint,60,.0001);
2938 if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2941 if (fDebugLevel>10){
2942 printf("A side fit parameters:\n");
2944 printf("\nC side fit parameters:\n");
2946 printf("\nAC side fit parameters:\n");
2953 //_____________________________________________________________________
2954 Double_t AliTPCCalibCE::SetBurstHnDrift()
2957 // Create a new THnSparse for the current burst
2958 // return the time of the current burst
2961 for(i=0; i<fTimeBursts.GetNrows(); ++i){
2962 if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2963 if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2964 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2965 return fTimeBursts(i);
2970 fArrHnDrift.AddAt(fHnDrift,i);
2971 fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
2976 //_____________________________________________________________________
2977 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
2980 // Write class to file
2981 // option can be specified in the dir option:
2983 // name=<objname>: the name of the calibration object in file will be <objname>
2984 // type=<type>: the saving type:
2985 // 0 - write the complte object
2986 // 1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
2987 // 2 - like 2, but in addition delete objects that will most probably not be used for calibration
2988 // 3 - store only calibration output, don't store the reference histograms
2989 // and THnSparse (requires Analyse called before)
2991 // NOTE: to read the object back, the ReadFromFile function should be used
2995 TString objName=GetName();
2999 TObjArray *arr=sDir.Tokenize(",");
3002 while ( (s=(TObjString*)next()) ){
3003 TString optString=s->GetString();
3004 optString.Remove(TString::kBoth,' ');
3005 if (optString.BeginsWith("name=")){
3006 objName=optString.ReplaceAll("name=","");
3008 if (optString.BeginsWith("type=")){
3009 optString.ReplaceAll("type=","");
3010 type=optString.Atoi();
3014 if (type==1||type==2) {
3015 //delete most probably not needed stuff
3016 //This requires Analyse to be called after reading the object from file
3017 fCalRocArrayT0.Delete();
3018 fCalRocArrayT0Err.Delete();
3019 fCalRocArrayQ.Delete();
3020 fCalRocArrayRMS.Delete();
3021 fCalRocArrayOutliers.Delete();
3023 if (type==2||type==3){
3024 fParamArrayEventPol1.Delete();
3025 fParamArrayEventPol2.Delete();
3028 TObjArray histoQArray(72);
3029 TObjArray histoT0Array(72);
3030 TObjArray histoRMSArray(72);
3031 TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3033 //save all large 2D histograms in separte pointers
3034 //to have a smaller memory print when saving the object
3035 if (type==1||type==2||type==3){
3036 for (Int_t i=0; i<72; ++i){
3037 histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3038 histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3039 histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3041 fHistoQArray.SetOwner(kFALSE);
3042 fHistoT0Array.SetOwner(kFALSE);
3043 fHistoRMSArray.SetOwner(kFALSE);
3044 fHistoQArray.Clear();
3045 fHistoT0Array.Clear();
3046 fHistoRMSArray.Clear();
3048 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3049 arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3051 fArrHnDrift.SetOwner(kFALSE);
3052 fArrHnDrift.Clear();
3056 TDirectory *backup = gDirectory;
3058 TFile f(filename,"recreate");
3059 this->Write(objName.Data());
3060 if (type==1||type==2) {
3061 histoQArray.Write("histoQArray",TObject::kSingleKey);
3062 histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3063 histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3064 arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3070 //move histograms back to the object
3071 if (type==1||type==2){
3072 for (Int_t i=0; i<72; ++i){
3073 fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3074 fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3075 fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3077 fHistoQArray.SetOwner(kTRUE);
3078 fHistoT0Array.SetOwner(kTRUE);
3079 fHistoRMSArray.SetOwner(kTRUE);
3081 for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3082 fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3084 fArrHnDrift.SetOwner(kTRUE);
3087 if ( backup ) backup->cd();
3089 //_____________________________________________________________________
3090 AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3093 // Read object from file
3094 // Handle properly if the histogram arrays were stored separately
3095 // call Analyse to make sure to have the calibration relevant information in the object
3099 if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3100 TList *l=f.GetListOfKeys();
3105 AliTPCCalibCE *ce=0x0;
3106 TObjArray *histoQArray=0x0;
3107 TObjArray *histoT0Array=0x0;
3108 TObjArray *histoRMSArray=0x0;
3109 TObjArray *arrHnDrift=0x0;
3111 while ( (key=(TKey*)next()) ){
3113 if ( o->IsA()==AliTPCCalibCE::Class() ){
3114 ce=(AliTPCCalibCE*)o;
3115 } else if ( o->IsA()==TObjArray::Class() ){
3116 TString name=key->GetName();
3117 if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3118 if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3119 if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3120 if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3125 //move histograms back to the object
3128 for (Int_t i=0; i<72; ++i){
3129 hist=(TH1*)histoQArray->UncheckedAt(i);
3130 if (hist) hist->SetDirectory(0x0);
3131 ce->fHistoQArray.AddAt(hist,i);
3133 ce->fHistoQArray.SetOwner(kTRUE);
3137 for (Int_t i=0; i<72; ++i){
3138 hist=(TH1*)histoT0Array->UncheckedAt(i);
3139 if (hist) hist->SetDirectory(0x0);
3140 ce->fHistoT0Array.AddAt(hist,i);
3142 ce->fHistoT0Array.SetOwner(kTRUE);
3146 for (Int_t i=0; i<72; ++i){
3147 hist=(TH1*)histoRMSArray->UncheckedAt(i);
3148 if (hist) hist->SetDirectory(0x0);
3149 ce->fHistoRMSArray.AddAt(hist,i);
3151 ce->fHistoRMSArray.SetOwner(kTRUE);
3155 for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3156 THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3157 ce->fArrHnDrift.AddAt(hSparse,i);