Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 :
17 : /// \class AliTPCCalROC
18 : ///
19 : /// Calibration base class for a single ROC
20 : /// Contains one float value per pad
21 : /// mapping of the pads taken form AliTPCROC
22 : ///
23 :
24 : // ROOT includes
25 : #include "TMath.h"
26 : #include "TClass.h"
27 : #include "TFile.h"
28 : #include "TH1F.h"
29 : #include "TH2F.h"
30 : #include "TArrayI.h"
31 :
32 : //
33 : //
34 : #include "AliTPCCalROC.h"
35 : #include "AliMathBase.h"
36 :
37 : #include "TRandom3.h" // only needed by the AliTPCCalROCTest() method
38 :
39 : /// \cond CLASSIMP
40 24 : ClassImp(AliTPCCalROC)
41 : /// \endcond
42 :
43 :
44 : //_____________________________________________________________________________
45 : AliTPCCalROC::AliTPCCalROC()
46 4896 : :TNamed(),
47 4896 : fSector(0),
48 4896 : fNChannels(0),
49 4896 : fNRows(0),
50 4896 : fkIndexes(0),
51 4896 : fData(0)
52 24480 : {
53 : //
54 : // Default constructor
55 : //
56 :
57 9792 : }
58 :
59 : //_____________________________________________________________________________
60 : AliTPCCalROC::AliTPCCalROC(UInt_t sector)
61 648 : :TNamed(),
62 648 : fSector(0),
63 648 : fNChannels(0),
64 648 : fNRows(0),
65 648 : fkIndexes(0),
66 648 : fData(0)
67 3240 : {
68 : /// Constructor that initializes a given sector
69 :
70 648 : fSector = sector;
71 1296 : fNChannels = AliTPCROC::Instance()->GetNChannels(fSector);
72 1296 : fNRows = AliTPCROC::Instance()->GetNRows(fSector);
73 1296 : fkIndexes = AliTPCROC::Instance()->GetRowIndexes(fSector);
74 1296 : fData = new Float_t[fNChannels];
75 648 : Reset();
76 :
77 1296 : }
78 :
79 : //_____________________________________________________________________________
80 : AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
81 0 : :TNamed(c),
82 0 : fSector(0),
83 0 : fNChannels(0),
84 0 : fNRows(0),
85 0 : fkIndexes(0),
86 0 : fData(0)
87 0 : {
88 : /// AliTPCCalROC copy constructor
89 :
90 0 : fSector = c.fSector;
91 0 : fNChannels = AliTPCROC::Instance()->GetNChannels(fSector);
92 0 : fNRows = AliTPCROC::Instance()->GetNRows(fSector);
93 0 : fkIndexes = AliTPCROC::Instance()->GetRowIndexes(fSector);
94 : //
95 0 : fData = new Float_t[fNChannels];
96 0 : for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
97 0 : }
98 : //____________________________________________________________________________
99 : AliTPCCalROC & AliTPCCalROC::operator =(const AliTPCCalROC & param)
100 : {
101 : /// assignment operator - dummy
102 :
103 0 : if (this == ¶m) return (*this);
104 0 : fSector = param.fSector;
105 0 : fNChannels = AliTPCROC::Instance()->GetNChannels(fSector);
106 0 : fNRows = AliTPCROC::Instance()->GetNRows(fSector);
107 0 : fkIndexes = AliTPCROC::Instance()->GetRowIndexes(fSector);
108 : //
109 0 : if (fData) delete [] fData;
110 0 : fData = new Float_t[fNChannels];
111 0 : for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata] = param.fData[idata];
112 0 : return (*this);
113 0 : }
114 :
115 :
116 : //_____________________________________________________________________________
117 : AliTPCCalROC::~AliTPCCalROC()
118 7344 : {
119 : /// AliTPCCalROC destructor
120 :
121 1224 : if (fData) {
122 2448 : delete [] fData;
123 1224 : fData = 0;
124 1224 : }
125 3672 : }
126 :
127 :
128 :
129 : void AliTPCCalROC::Streamer(TBuffer &R__b)
130 : {
131 : /// Stream an object of class AliTPCCalROC.
132 :
133 14688 : if (R__b.IsReading()) {
134 9792 : AliTPCCalROC::Class()->ReadBuffer(R__b, this);
135 4896 : fkIndexes = AliTPCROC::Instance()->GetRowIndexes(fSector);
136 4896 : } else {
137 0 : AliTPCCalROC::Class()->WriteBuffer(R__b,this);
138 : }
139 4896 : }
140 :
141 : void AliTPCCalROC::Print(Option_t* /*option*/) const{
142 : //
143 : // Print summary content
144 : //
145 0 : printf("|ROC%d|\t%f|\t%f|\t%f|\n",fSector, GetMean(),GetMedian(),GetRMS());
146 0 : }
147 :
148 :
149 :
150 : //
151 :
152 : Bool_t AliTPCCalROC::MedianFilter(Int_t deltaRow, Int_t deltaPad, AliTPCCalROC* outlierROC, Bool_t doEdge){
153 : /// Modify content of the object - raplace value by median in neighorhood
154 :
155 0 : Float_t *newBuffer=new Float_t[fNChannels] ;
156 0 : Double_t *cacheBuffer=new Double_t[fNChannels];
157 : //
158 0 : for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
159 0 : Int_t nPads=GetNPads(iRow); // number of rows in current row
160 0 : for (Int_t iPad=0; iPad<nPads; iPad++){
161 0 : Double_t value=GetValue(iRow,iPad);
162 : Int_t counter=0;
163 : //
164 0 : for (Int_t dRow=-deltaRow; dRow<=deltaRow; dRow++){
165 0 : Int_t jRow=iRow+dRow; //take the row - mirror if ouside of range
166 : Float_t sign0=1.;
167 0 : if (jRow<0) sign0=-1.;
168 0 : if (UInt_t(jRow)>=fNRows) sign0=-1.;
169 0 : Int_t jPads= GetNPads(iRow+sign0*dRow);
170 0 : Int_t offset=(nPads-jPads)/2.;
171 : //
172 0 : for (Int_t dPad=-deltaPad; dPad<=deltaPad; dPad++){
173 : Float_t sign=sign0;
174 0 : Int_t jPad=iPad+dPad+offset; //take the pad - mirror if ouside of range
175 : Int_t kRow=jRow;
176 0 : if (jPad<0) sign=-1;
177 0 : if (jPad>=jPads) sign=-1;
178 0 : if (sign<0){
179 0 : kRow=iRow-dRow;
180 0 : jPad=iPad-dPad+offset;
181 0 : if (!doEdge) continue;
182 : }
183 0 : if (IsInRange(UInt_t(kRow),jPad)){
184 0 : Bool_t isOutlier=(outlierROC==NULL)?kFALSE:outlierROC->GetValue(kRow,jPad)>0;
185 0 : if (!isOutlier){
186 0 : cacheBuffer[counter]=sign*(GetValue(kRow,jPad)-value);
187 0 : counter++;
188 0 : }
189 0 : }
190 0 : }
191 : }
192 0 : newBuffer[fkIndexes[iRow]+iPad]=0.;
193 0 : if (counter>1) newBuffer[fkIndexes[iRow]+iPad] = TMath::Median(counter, cacheBuffer)+value;
194 : }
195 : }
196 0 : memcpy(fData, newBuffer,GetNchannels()*sizeof(Float_t));
197 0 : delete []newBuffer;
198 0 : delete []cacheBuffer;
199 0 : return kTRUE;
200 : }
201 :
202 :
203 :
204 : Bool_t AliTPCCalROC::LTMFilter(Int_t deltaRow, Int_t deltaPad, Float_t fraction, Int_t type, AliTPCCalROC* outlierROC, Bool_t doEdge){
205 : /// Modify content of the class
206 : /// write LTM mean or median
207 :
208 0 : if (fraction<0 || fraction>1) return kFALSE;
209 0 : Float_t *newBuffer=new Float_t[fNChannels] ;
210 0 : Double_t *cacheBuffer=new Double_t[fNChannels];
211 : //
212 0 : for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
213 0 : Int_t nPads=GetNPads(iRow); // number of rows in current row
214 0 : for (Int_t iPad=0; iPad<nPads; iPad++){
215 0 : Double_t value=GetValue(iRow,iPad);
216 : Int_t counter=0;
217 : //
218 0 : for (Int_t dRow=-deltaRow; dRow<=deltaRow; dRow++){
219 0 : Int_t jRow=iRow+dRow; //take the row - mirror if ouside of range
220 : Float_t sign0=1.;
221 0 : if (jRow<0) sign0=-1.;
222 0 : if (UInt_t(jRow)>=fNRows) sign0=-1.;
223 0 : Int_t jPads= GetNPads(iRow+sign0*dRow);
224 0 : Int_t offset=(nPads-jPads)/2.;
225 : //
226 0 : for (Int_t dPad=-deltaPad; dPad<=deltaPad; dPad++){
227 : Float_t sign=sign0;
228 0 : Int_t jPad=iPad+dPad+offset; //take the pad - mirror if ouside of range
229 : Int_t kRow=jRow;
230 0 : if (jPad<0) sign=-1;
231 0 : if (jPad>=jPads) sign=-1;
232 0 : if (sign<0){
233 0 : if (!doEdge) continue;
234 0 : kRow=iRow-dRow;
235 0 : jPad=iPad-dPad+offset;
236 0 : }
237 0 : if (IsInRange(UInt_t(kRow),jPad)){
238 0 : Bool_t isOutlier=(outlierROC==NULL)?kFALSE:outlierROC->GetValue(kRow,jPad)>0;
239 0 : if (!isOutlier){
240 0 : cacheBuffer[counter]=sign*(GetValue(kRow,jPad)-value);
241 0 : counter++;
242 0 : }
243 0 : }
244 0 : }
245 : }
246 0 : Double_t mean=0,rms=0;
247 0 : if (TMath::Nint(fraction*Double_t(counter))>1 ) AliMathBase::EvaluateUni(counter,cacheBuffer,mean,rms,TMath::Nint(fraction*Double_t(counter)));
248 0 : mean+=value;
249 0 : newBuffer[fkIndexes[iRow]+iPad] = (type==0)? mean:rms;
250 0 : }
251 : }
252 0 : memcpy(fData, newBuffer,GetNchannels()*sizeof(Float_t));
253 0 : delete []newBuffer;
254 0 : delete []cacheBuffer;
255 : return kTRUE;
256 0 : }
257 :
258 : Bool_t AliTPCCalROC::Convolute(Double_t sigmaPad, Double_t sigmaRow, AliTPCCalROC*outlierROC, TF1 */*fpad*/, TF1 */*frow*/){
259 : /// convolute the calibration with function fpad,frow
260 : /// in range +-4 sigma
261 :
262 0 : Float_t *newBuffer=new Float_t[fNChannels] ;
263 : //
264 0 : for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
265 0 : Int_t nPads=GetNPads(iRow); // number of rows in current row
266 0 : for (Int_t iPad=0; iPad<nPads; iPad++){
267 0 : Int_t jRow0=TMath::Max(TMath::Nint(iRow-sigmaRow*4.),0);
268 0 : Int_t jRow1=TMath::Min(TMath::Nint(iRow+sigmaRow*4.),Int_t(fNRows));
269 0 : Int_t jPad0=TMath::Max(TMath::Nint(iPad-sigmaPad*4.),0);
270 0 : Int_t jPad1=TMath::Min(TMath::Nint(iRow+sigmaPad*4.),Int_t(nPads));
271 : //
272 : Double_t sumW=0;
273 : Double_t sumCal=0;
274 0 : for (Int_t jRow=jRow0; jRow<=jRow1; jRow++){
275 0 : for (Int_t jPad=jPad0; jPad<=jPad1; jPad++){
276 0 : if (!IsInRange(jRow,jPad)) continue;
277 0 : Bool_t isOutlier=(outlierROC==NULL)?kFALSE:outlierROC->GetValue(jRow,jPad)>0;
278 0 : if (!isOutlier){
279 0 : Double_t weight= TMath::Gaus(jPad,iPad,sigmaPad)*TMath::Gaus(jRow,iRow,sigmaRow);
280 0 : sumCal+=weight*GetValue(jRow,jPad);
281 0 : sumW+=weight;
282 0 : }
283 0 : }
284 : }
285 0 : if (sumW>0){
286 0 : sumCal/=sumW;
287 0 : newBuffer[fkIndexes[iRow]+iPad]=sumCal;
288 0 : }
289 : }
290 : }
291 0 : memcpy(fData, newBuffer,GetNchannels()*sizeof(Float_t));
292 0 : delete []newBuffer;
293 0 : return kTRUE;
294 : }
295 :
296 :
297 : //
298 :
299 :
300 :
301 :
302 :
303 : // algebra fuctions:
304 :
305 : void AliTPCCalROC::Add(Float_t c1){
306 : /// add c1 to each channel of the ROC
307 :
308 0 : for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
309 0 : }
310 :
311 :
312 : void AliTPCCalROC::Multiply(Float_t c1){
313 : /// multiply each channel of the ROC with c1
314 :
315 0 : for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
316 0 : }
317 :
318 :
319 : void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
320 : /// multiply AliTPCCalROC roc by c1 and add each channel to the coresponing channel in the ROC
321 : /// - pad by pad
322 :
323 0 : if (!roc) return;
324 0 : for (UInt_t idata = 0; idata< fNChannels; idata++){
325 0 : fData[idata]+=roc->fData[idata]*c1;
326 : }
327 0 : }
328 :
329 :
330 : void AliTPCCalROC::Multiply(const AliTPCCalROC* roc) {
331 : /// multiply each channel of the ROC with the coresponding channel of 'roc'
332 : /// - pad by pad -
333 :
334 0 : if (!roc) return;
335 0 : for (UInt_t idata = 0; idata< fNChannels; idata++){
336 0 : fData[idata]*=roc->fData[idata];
337 : }
338 0 : }
339 :
340 :
341 : void AliTPCCalROC::Divide(const AliTPCCalROC* roc) {
342 : /// divide each channel of the ROC by the coresponding channel of 'roc'
343 : /// - pad by pad -
344 :
345 0 : if (!roc) return;
346 : Float_t kEpsilon=0.00000000000000001;
347 0 : for (UInt_t idata = 0; idata< fNChannels; idata++){
348 0 : if (TMath::Abs(roc->fData[idata])>kEpsilon)
349 0 : fData[idata]/=roc->fData[idata];
350 : }
351 0 : }
352 :
353 : void AliTPCCalROC::Reset()
354 : {
355 : /// reset to ZERO
356 :
357 1296 : memset(fData,0,sizeof(Float_t)*fNChannels); // set all values to 0
358 648 : }
359 :
360 : //____________________________________________________________
361 : Double_t AliTPCCalROC::GetStats(EStatType statType, AliTPCCalROC *const outlierROC, EPadType padType) const
362 : {
363 : /// returns the mean or median or RMS value depending on the statType
364 : /// of the ROC or medium or long pad region
365 : /// pads with value != 0 in outlierROC are not used for the calculation
366 : /// return 0 if no data is accepted by the outlier cuts
367 :
368 0 : if(fSector<36) padType=kAll;
369 0 : if (!outlierROC) {
370 0 : if(padType==kAll) {
371 0 : if (statType==kMean) return TMath::Mean (fNChannels, fData);
372 0 : else if(statType==kMedian) return TMath::Median (fNChannels, fData);
373 0 : else if(statType==kRMS) return TMath::RMS (fNChannels, fData);
374 0 : else if(statType==kMinElement) return TMath::MinElement(fNChannels, fData);
375 0 : else if(statType==kMaxElement) return TMath::MaxElement(fNChannels, fData);
376 : }
377 0 : else if(padType==kOROCmedium) {
378 0 : if (statType==kMean) return TMath::Mean (fkIndexes[64], fData);
379 0 : else if(statType==kMedian) return TMath::Median (fkIndexes[64], fData);
380 0 : else if(statType==kRMS) return TMath::RMS (fkIndexes[64], fData);
381 0 : else if(statType==kMinElement) return TMath::MinElement(fkIndexes[64], fData);
382 0 : else if(statType==kMaxElement) return TMath::MaxElement(fkIndexes[64], fData);
383 : }
384 0 : else if(padType==kOROClong) {
385 0 : const Int_t offset=fkIndexes[64];
386 0 : if (statType==kMean) return TMath::Mean (fNChannels-offset, fData+offset);
387 0 : else if(statType==kMedian) return TMath::Median (fNChannels-offset, fData+offset);
388 0 : else if(statType==kRMS) return TMath::RMS (fNChannels-offset, fData+offset);
389 0 : else if(statType==kMinElement) return TMath::MinElement(fNChannels-offset, fData+offset);
390 0 : else if(statType==kMaxElement) return TMath::MaxElement(fNChannels-offset, fData+offset);
391 0 : }
392 : }//end of no outlierROC
393 :
394 0 : Float_t ddata[fNChannels];
395 :
396 : Int_t nPoints = 0;
397 : UInt_t indexMin = 0;
398 0 : UInt_t indexMax = fNChannels;
399 0 : if(padType==kOROCmedium) indexMax=fkIndexes[64];
400 0 : else if(padType==kOROClong) indexMin=fkIndexes[64];
401 :
402 0 : for (UInt_t i=indexMin;i<indexMax;i++) {
403 0 : if (outlierROC->GetValue(i)>1e-20) continue;
404 0 : ddata[nPoints]= fData[i];
405 0 : nPoints++;
406 0 : }
407 : Double_t value = 0;
408 0 : if(nPoints>0){
409 0 : if (statType==kMean) value = TMath::Mean (nPoints, ddata);
410 0 : else if(statType==kMedian) value = TMath::Median (nPoints, ddata);
411 0 : else if(statType==kRMS) value = TMath::RMS (nPoints, ddata);
412 0 : else if(statType==kMinElement) value = TMath::MinElement(nPoints, ddata);
413 0 : else if(statType==kMaxElement) value = TMath::MaxElement(nPoints, ddata);
414 : }
415 : return value;
416 0 : }
417 :
418 : //_____________________________________________________________________________
419 : Double_t AliTPCCalROC::GetMean(AliTPCCalROC *const outlierROC,EPadType padType) const
420 : {
421 : /// returns the mean value of the ROC
422 : /// pads with value != 0 in outlierROC are not used for the calculation
423 : /// return 0 if no data is accepted by the outlier cuts
424 :
425 0 : return GetStats(kMean,outlierROC,padType);
426 : }
427 :
428 : //_____________________________________________________________________________
429 : Double_t AliTPCCalROC::GetMedian(AliTPCCalROC *const outlierROC,EPadType padType) const
430 : {
431 : /// returns the median value of the ROC
432 : /// pads with value != 0 in outlierROC are not used for the calculation
433 : /// return 0 if no data is accepted by the outlier cuts
434 :
435 0 : return GetStats(kMedian,outlierROC,padType);
436 : }
437 :
438 : //_____________________________________________________________________________
439 : Double_t AliTPCCalROC::GetRMS(AliTPCCalROC *const outlierROC,EPadType padType) const
440 : {
441 : /// returns the RMS value of the ROC
442 : /// pads with value != 0 in outlierROC are not used for the calculation
443 : /// return 0 if no data is accepted by the outlier cuts
444 :
445 0 : return GetStats(kRMS,outlierROC,padType);
446 : }
447 :
448 : //_____________________________________________________________________________
449 : Double_t AliTPCCalROC::GetMinElement(AliTPCCalROC *const outlierROC,EPadType padType) const
450 : {
451 : /// returns the MinElement value of the ROC
452 : /// pads with value != 0 in outlierROC are not used for the calculation
453 : /// return 0 if no data is accepted by the outlier cuts
454 :
455 0 : return GetStats(kMinElement,outlierROC,padType);
456 : }
457 :
458 : //_____________________________________________________________________________
459 : Double_t AliTPCCalROC::GetMaxElement(AliTPCCalROC *const outlierROC,EPadType padType) const
460 : {
461 : /// returns the MaxElement value of the ROC
462 : /// pads with value != 0 in outlierROC are not used for the calculation
463 : /// return 0 if no data is accepted by the outlier cuts
464 :
465 0 : return GetStats(kMaxElement,outlierROC,padType);
466 : }
467 :
468 : //_____________________________________________________________________________
469 : Double_t AliTPCCalROC::GetLTM(Double_t *const sigma, Double_t fraction, AliTPCCalROC *const outlierROC, EPadType padType){
470 : /// returns the LTM and sigma
471 : /// pads with value != 0 in outlierROC are not used for the calculation
472 : /// return 0 if no data is accepted by the outlier cuts
473 : /// LTM for different padType
474 :
475 0 : if(fSector<36) padType=kAll;
476 0 : Double_t ddata[fNChannels];
477 : UInt_t nPoints = 0;
478 : UInt_t indexMin = 0;
479 0 : UInt_t indexMax = fNChannels;
480 0 : if(padType==kOROCmedium) indexMax=fkIndexes[64];
481 0 : else if(padType==kOROClong) indexMin=fkIndexes[64];
482 0 : for (UInt_t i=indexMin;i<indexMax;i++) {
483 0 : if (outlierROC && (outlierROC->GetValue(i) >1e-20)) continue;
484 0 : ddata[nPoints]= fData[i];
485 0 : nPoints++;
486 0 : }
487 :
488 0 : Double_t ltm =0, lsigma=0;
489 0 : if(nPoints>0) {
490 0 : Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
491 0 : AliMathBase::EvaluateUni(nPoints,ddata, ltm, lsigma, hh);
492 0 : if (sigma) *sigma=lsigma;
493 0 : }
494 :
495 0 : return ltm;
496 0 : }
497 :
498 : //_____________________________________________________________________________
499 : TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
500 : /// make 1D histo
501 : /// type -1 = user defined range
502 : /// 0 = nsigma cut nsigma=min
503 : /// 1 = delta cut around median delta=min
504 :
505 0 : if (type>=0){
506 0 : if (type==0){
507 : // nsigma range
508 0 : Float_t mean = GetMean();
509 0 : Float_t sigma = GetRMS();
510 0 : Float_t nsigma = TMath::Abs(min);
511 0 : min = mean-nsigma*sigma;
512 0 : max = mean+nsigma*sigma;
513 0 : }
514 0 : if (type==1){
515 : // fixed range
516 0 : Float_t mean = GetMedian();
517 : Float_t delta = min;
518 0 : min = mean-delta;
519 0 : max = mean+delta;
520 0 : }
521 0 : if (type==2){
522 : //
523 : // LTM mean +- nsigma
524 : //
525 0 : Double_t sigma;
526 0 : Float_t mean = GetLTM(&sigma,max);
527 0 : sigma*=min;
528 0 : min = mean-sigma;
529 0 : max = mean+sigma;
530 0 : }
531 : }
532 0 : TString name=Form("%s ROC 1D%d",GetTitle(),fSector);
533 0 : TH1F * his = new TH1F(name.Data(),name.Data(),100, min,max);
534 0 : for (UInt_t irow=0; irow<fNRows; irow++){
535 0 : UInt_t npads = (Int_t)GetNPads(irow);
536 0 : for (UInt_t ipad=0; ipad<npads; ipad++){
537 0 : his->Fill(GetValue(irow,ipad));
538 : }
539 : }
540 : return his;
541 0 : }
542 :
543 :
544 :
545 : TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
546 : /// make 2D histo
547 : /// type -1 = user defined range
548 : /// 0 = nsigma cut nsigma=min
549 : /// 1 = delta cut around median delta=min
550 :
551 0 : if (type>=0){
552 0 : if (type==0){
553 : // nsigma range
554 0 : Float_t mean = GetMean();
555 0 : Float_t sigma = GetRMS();
556 0 : Float_t nsigma = TMath::Abs(min);
557 0 : min = mean-nsigma*sigma;
558 0 : max = mean+nsigma*sigma;
559 0 : }
560 0 : if (type==1){
561 : // fixed range
562 0 : Float_t mean = GetMedian();
563 : Float_t delta = min;
564 0 : min = mean-delta;
565 0 : max = mean+delta;
566 0 : }
567 0 : if (type==2){
568 0 : Double_t sigma;
569 0 : Float_t mean = GetLTM(&sigma,max);
570 0 : sigma*=min;
571 0 : min = mean-sigma;
572 0 : max = mean+sigma;
573 :
574 0 : }
575 : }
576 : UInt_t maxPad = 0;
577 0 : for (UInt_t irow=0; irow<fNRows; irow++){
578 0 : if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
579 : }
580 :
581 0 : TString name=Form("%s ROC%d",GetTitle(),fSector);
582 0 : TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
583 0 : for (UInt_t irow=0; irow<fNRows; irow++){
584 0 : UInt_t npads = (Int_t)GetNPads(irow);
585 0 : for (UInt_t ipad=0; ipad<npads; ipad++){
586 0 : his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
587 : }
588 : }
589 0 : his->SetMaximum(max);
590 0 : his->SetMinimum(min);
591 : return his;
592 0 : }
593 :
594 : TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
595 : /// Make Histogram with outliers
596 : /// mode = 0 - sigma cut used
597 : /// mode = 1 - absolute cut used
598 : /// fraction - fraction of values used to define sigma
599 : /// delta - in mode 0 - nsigma cut
600 : /// mode 1 - delta cut
601 :
602 0 : Double_t sigma;
603 0 : Float_t mean = GetLTM(&sigma,fraction);
604 0 : if (type==0) delta*=sigma;
605 : UInt_t maxPad = 0;
606 0 : for (UInt_t irow=0; irow<fNRows; irow++){
607 0 : if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
608 : }
609 :
610 0 : TString name=Form("%s ROC Outliers%d",GetTitle(),fSector);
611 0 : TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
612 0 : for (UInt_t irow=0; irow<fNRows; irow++){
613 0 : UInt_t npads = (Int_t)GetNPads(irow);
614 0 : for (UInt_t ipad=0; ipad<npads; ipad++){
615 0 : if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
616 0 : his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
617 : }
618 : }
619 : return his;
620 0 : }
621 :
622 :
623 :
624 : void AliTPCCalROC::Draw(Option_t* opt){
625 : /// create histogram with values and draw it
626 :
627 : TH1 * his=0;
628 0 : TString option=opt;
629 0 : option.ToUpper();
630 0 : if (option.Contains("1D")){
631 0 : his = MakeHisto1D();
632 0 : }
633 : else{
634 0 : his = MakeHisto2D();
635 : }
636 0 : his->Draw(option);
637 0 : }
638 :
639 :
640 :
641 :
642 :
643 : void AliTPCCalROC::Test() {
644 : /// example function to show functionality and test AliTPCCalROC
645 :
646 : Float_t kEpsilon=0.00001;
647 :
648 : // create CalROC with defined values
649 0 : AliTPCCalROC roc0(0);
650 0 : for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
651 0 : for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
652 0 : Float_t value = irow+ipad/1000.;
653 0 : roc0.SetValue(irow,ipad,value);
654 : }
655 : }
656 :
657 : // copy CalROC, readout values and compare to original
658 0 : AliTPCCalROC roc1(roc0);
659 0 : for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
660 0 : for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
661 0 : Float_t value = irow+ipad/1000.;
662 0 : if (roc1.GetValue(irow,ipad)!=value){
663 0 : printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
664 : }
665 : }
666 : }
667 :
668 : // write original CalROC to file
669 0 : TFile f("calcTest.root","recreate");
670 0 : roc0.Write("Roc0");
671 0 : AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
672 0 : f.Close();
673 :
674 : // read CalROC from file and compare to original CalROC
675 0 : for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
676 0 : if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
677 0 : printf("NPads - Read/Write error\trow=%d\n",irow);
678 0 : for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
679 0 : Float_t value = irow+ipad/1000.;
680 0 : if (roc2->GetValue(irow,ipad)!=value){
681 0 : printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
682 : }
683 : }
684 : }
685 :
686 : //
687 : // Algebra Tests
688 : //
689 :
690 : // Add constant
691 0 : AliTPCCalROC roc3(roc0);
692 0 : roc3.Add(1.5);
693 0 : for (UInt_t irow = 0; irow <roc3.GetNrows(); irow++){
694 0 : for (UInt_t ipad = 0; ipad <roc3.GetNPads(irow); ipad++){
695 0 : Float_t value = irow+ipad/1000. + 1.5;
696 0 : if (TMath::Abs(roc3.GetValue(irow,ipad)-value) > kEpsilon){
697 0 : printf("Add constant - error\trow=%d\tpad=%d\n",irow,ipad);
698 : }
699 : }
700 : }
701 :
702 : // Add another CalROC
703 0 : AliTPCCalROC roc4(roc0);
704 0 : roc4.Add(&roc0, -1.5);
705 0 : for (UInt_t irow = 0; irow <roc4.GetNrows(); irow++){
706 0 : for (UInt_t ipad = 0; ipad <roc4.GetNPads(irow); ipad++){
707 0 : Float_t value = irow+ipad/1000. - 1.5 * (irow+ipad/1000.);
708 0 : if (TMath::Abs(roc4.GetValue(irow,ipad)-value) > kEpsilon){
709 0 : printf("Add CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
710 : }
711 : }
712 : }
713 :
714 : // Multiply with constant
715 0 : AliTPCCalROC roc5(roc0);
716 0 : roc5.Multiply(-1.4);
717 0 : for (UInt_t irow = 0; irow <roc5.GetNrows(); irow++){
718 0 : for (UInt_t ipad = 0; ipad <roc5.GetNPads(irow); ipad++){
719 0 : Float_t value = (irow+ipad/1000.) * (-1.4);
720 0 : if (TMath::Abs(roc5.GetValue(irow,ipad)-value) > kEpsilon){
721 0 : printf("Multiply with constant - error\trow=%d\tpad=%d\n",irow,ipad);
722 : }
723 : }
724 : }
725 :
726 : // Multiply another CalROC
727 0 : AliTPCCalROC roc6(roc0);
728 0 : roc6.Multiply(&roc0);
729 0 : for (UInt_t irow = 0; irow <roc6.GetNrows(); irow++){
730 0 : for (UInt_t ipad = 0; ipad <roc6.GetNPads(irow); ipad++){
731 0 : Float_t value = (irow+ipad/1000.) * (irow+ipad/1000.);
732 0 : if (TMath::Abs(roc6.GetValue(irow,ipad)-value) > kEpsilon*100){
733 0 : printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
734 : }
735 : }
736 : }
737 :
738 :
739 : // Divide by CalROC
740 0 : AliTPCCalROC roc7(roc0);
741 0 : roc7.Divide(&roc0);
742 0 : for (UInt_t irow = 0; irow <roc7.GetNrows(); irow++){
743 0 : for (UInt_t ipad = 0; ipad <roc7.GetNPads(irow); ipad++){
744 : Float_t value = 1.;
745 0 : if (irow+ipad == 0) value = roc0.GetValue(irow,ipad);
746 0 : if (TMath::Abs(roc7.GetValue(irow,ipad)-value) > kEpsilon){
747 0 : printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
748 : }
749 : }
750 : }
751 :
752 : //
753 : // Statistics Test
754 : //
755 :
756 : // create CalROC with defined values
757 0 : TRandom3 rnd(0);
758 0 : AliTPCCalROC sroc0(0);
759 0 : for (UInt_t ichannel = 0; ichannel < sroc0.GetNchannels(); ichannel++){
760 0 : Float_t value = rnd.Gaus(10., 2.);
761 0 : sroc0.SetValue(ichannel,value);
762 : }
763 :
764 0 : printf("Mean (should be close to 10): %f\n", sroc0.GetMean());
765 0 : printf("RMS (should be close to 2): %f\n", sroc0.GetRMS());
766 0 : printf("Median (should be close to 10): %f\n", sroc0.GetMedian());
767 0 : printf("LTM (should be close to 10): %f\n", sroc0.GetLTM());
768 :
769 : //AliTPCCalROC* sroc1 = sroc0.LocalFit(4, 8);
770 :
771 : //delete sroc1;
772 :
773 : // std::cout << TMath::Abs(roc5.GetValue(irow,ipad)-value) << std::endl;
774 0 : }
775 :
776 :
777 : AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
778 : /// MakeLocalFit - smoothing
779 : /// returns a AliTPCCalROC with smoothed data
780 : /// rowRadius and padRadius specifies a window around a given pad.
781 : /// The data of this window are fitted with a parabolic function.
782 : /// This function is evaluated at the pad's position.
783 : /// At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
784 : /// rowRadius - radius - rows to be used for smoothing
785 : /// padradius - radius - pads to be used for smoothing
786 : /// ROCoutlier - map of outliers - pads not to be used for local smoothing
787 : /// robust - robust method of fitting - (much slower)
788 : /// chi2Threshold: Threshold for chi2 when EvalRobust is called
789 : /// robustFraction: Fraction of data that will be used in EvalRobust
790 :
791 0 : AliTPCCalROC * xROCfitted = new AliTPCCalROC(fSector);
792 0 : TLinearFitter fitterQ(6,"hyp5");
793 : // TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
794 0 : fitterQ.StoreData(kTRUE);
795 0 : for (UInt_t row=0; row < GetNrows(); row++) {
796 : //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
797 0 : for (UInt_t pad=0; pad < GetNPads(row); pad++)
798 0 : xROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust, chi2Threshold, robustFraction));
799 : }
800 : return xROCfitted;
801 0 : }
802 :
803 :
804 : Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC *const ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
805 : /// AliTPCCalROC::GetNeighbourhoodValue - smoothing - PRIVATE
806 : /// in this function the fit for LocalFit is done
807 :
808 0 : fitterQ->ClearPoints();
809 0 : TVectorD fitParam(6);
810 : Int_t npoints=0;
811 0 : Double_t xx[6];
812 : Float_t dlx, dly;
813 0 : Float_t lPad[3] = {0};
814 0 : Float_t localXY[3] = {0};
815 :
816 0 : AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
817 0 : tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad); // calculate position lPad by pad and row number
818 :
819 0 : TArrayI *neighbourhoodRows = 0;
820 0 : TArrayI *neighbourhoodPads = 0;
821 :
822 : //std::cerr << "Trying to get neighbourhood for row " << row << ", pad " << pad << std::endl;
823 0 : GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
824 : //std::cerr << "Got neighbourhood for row " << row << ", pad " << pad << std::endl;
825 :
826 : Int_t r, p;
827 0 : for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
828 0 : r = neighbourhoodRows->At(i);
829 0 : p = neighbourhoodPads->At(i);
830 0 : if (r == -1 || p == -1) continue; // window is bigger than ROC
831 0 : tpcROCinstance->GetPositionLocal(fSector, r, p, localXY); // calculate position localXY by pad and row number
832 0 : dlx = lPad[0] - localXY[0];
833 0 : dly = lPad[1] - localXY[1];
834 : //xx[0] = 1;
835 0 : xx[1] = dlx;
836 0 : xx[2] = dly;
837 0 : xx[3] = dlx*dlx;
838 0 : xx[4] = dly*dly;
839 0 : xx[5] = dlx*dly;
840 0 : if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) {
841 0 : fitterQ->AddPoint(&xx[1], GetValue(r, p), 1);
842 0 : npoints++;
843 0 : }
844 : }
845 :
846 0 : delete neighbourhoodRows;
847 0 : delete neighbourhoodPads;
848 :
849 0 : if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
850 : // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
851 0 : return 0.; // for diagnostic
852 : }
853 0 : fitterQ->Eval();
854 0 : fitterQ->GetParameters(fitParam);
855 : Float_t chi2Q = 0;
856 0 : if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
857 : //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
858 0 : if (robust && chi2Q > chi2Threshold) {
859 : //std::cout << "robust fitter called... " << std::endl;
860 0 : fitterQ->EvalRobust(robustFraction);
861 0 : fitterQ->GetParameters(fitParam);
862 : }
863 0 : Double_t value = fitParam[0];
864 :
865 : //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q << std::endl;
866 : return value;
867 0 : }
868 :
869 :
870 :
871 :
872 : void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
873 : /// AliTPCCalROC::GetNeighbourhood - PRIVATE
874 : /// in this function the window for LocalFit is determined
875 :
876 0 : rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
877 0 : padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
878 :
879 0 : Int_t rmin = row - rRadius;
880 0 : UInt_t rmax = row + rRadius;
881 :
882 : // if window goes out of ROC
883 0 : if (rmin < 0) {
884 0 : rmax = rmax - rmin;
885 : rmin = 0;
886 0 : }
887 0 : if (rmax >= GetNrows()) {
888 0 : rmin = rmin - (rmax - GetNrows()+1);
889 0 : rmax = GetNrows() - 1;
890 0 : if (rmin < 0 ) rmin = 0; // if the window is bigger than the ROC
891 : }
892 :
893 : Int_t pmin, pmax;
894 : Int_t i = 0;
895 :
896 0 : for (UInt_t r = rmin; r <= rmax; r++) {
897 0 : pmin = pad - pRadius;
898 0 : pmax = pad + pRadius;
899 0 : if (pmin < 0) {
900 0 : pmax = pmax - pmin;
901 : pmin = 0;
902 0 : }
903 0 : if (pmax >= (Int_t)GetNPads(r)) {
904 0 : pmin = pmin - (pmax - GetNPads(r)+1);
905 0 : pmax = GetNPads(r) - 1;
906 0 : if (pmin < 0 ) pmin = 0; // if the window is bigger than the ROC
907 : }
908 0 : for (Int_t p = pmin; p <= pmax; p++) {
909 0 : (*rowArray)[i] = r;
910 0 : (*padArray)[i] = p;
911 0 : i++;
912 : }
913 : }
914 0 : for (Int_t j = i; j < rowArray->GetSize(); j++){ // unused padArray-entries, in the case that the window is bigger than the ROC
915 : //std::cout << "trying to write -1" << std::endl;
916 0 : (*rowArray)[j] = -1;
917 0 : (*padArray)[j] = -1;
918 : //std::cout << "writing -1" << std::endl;
919 : }
920 0 : }
921 :
922 :
923 :
924 : void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t & chi2, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction, Double_t err, EPadType padType){
925 : /// Makes a GlobalFit for the given sector and return fit-parameters, covariance and chi2
926 : /// update of GlobalFit in 2015 to have different gain in OROC medium and long pads
927 : /// The origin of the fit function is the center of the ROC!
928 : /// fitType == 0: fit plane function
929 : /// fitType == 1: fit parabolic function
930 : /// ROCoutliers - pads with value !=0 are not used in fitting procedure
931 : /// chi2Threshold: Threshold for chi2 when EvalRobust is called
932 : /// robustFraction: Fraction of data that will be used in EvalRobust
933 : /// err: error of the data points
934 :
935 0 : if(fSector<36) padType=kAll;
936 :
937 : TLinearFitter* fitterG = 0;
938 0 : Double_t xx[6];
939 :
940 0 : if (fitType == 1) {
941 0 : fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
942 0 : fitParam.ResizeTo(6);
943 0 : covMatrix.ResizeTo(6,6);
944 0 : } else {
945 0 : fitterG = new TLinearFitter(3,"x0++x1++x2");
946 0 : fitParam.ResizeTo(3);
947 0 : covMatrix.ResizeTo(3,3);
948 : }
949 0 : fitterG->StoreData(kTRUE);
950 0 : fitterG->ClearPoints();
951 : Int_t npoints=0;
952 :
953 : Float_t dlx, dly;
954 0 : Float_t centerPad[3] = {0};
955 0 : Float_t localXY[3] = {0};
956 :
957 0 : AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
958 0 : if(padType==kAll)
959 0 : tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad); // calculate center of ROC
960 0 : else if(padType==kOROCmedium)
961 0 : tpcROCinstance->GetPositionLocal(fSector, 64/2, GetNPads(64/2)/2, centerPad); // calculate center of ROC medium pads
962 0 : else if(padType==kOROClong)
963 0 : tpcROCinstance->GetPositionLocal(fSector, (GetNrows()-64)/2+64, GetNPads((GetNrows()-64)/2+64)/2, centerPad); // calculate center of ROC long pads
964 :
965 : UInt_t irowMin = 0;
966 0 : UInt_t irowMax = GetNrows();
967 0 : if(padType==kOROCmedium) irowMax = 64;
968 0 : else if(padType==kOROClong) irowMin = 64;
969 :
970 : // loop over all channels and read data into fitterG
971 0 : for (UInt_t irow = irowMin; irow < irowMax; irow++) {
972 0 : for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
973 : // fill fitterG
974 0 : if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 0) continue;
975 0 : tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY); // calculate position localXY by pad and row number
976 0 : dlx = localXY[0] - centerPad[0];
977 0 : dly = localXY[1] - centerPad[1];
978 0 : xx[0] = 1;
979 0 : xx[1] = dlx;
980 0 : xx[2] = dly;
981 0 : xx[3] = dlx*dlx;
982 0 : xx[4] = dly*dly;
983 0 : xx[5] = dlx*dly;
984 0 : npoints++;
985 0 : fitterG->AddPoint(xx, GetValue(irow, ipad), err);
986 0 : }
987 : }
988 0 : if(npoints>10) { // make sure there is something to fit
989 0 : fitterG->Eval();
990 0 : fitterG->GetParameters(fitParam);
991 0 : fitterG->GetCovarianceMatrix(covMatrix);
992 0 : if (fitType == 1)
993 0 : chi2 = fitterG->GetChisquare()/(npoints-6.);
994 0 : else chi2 = fitterG->GetChisquare()/(npoints-3.);
995 0 : if (robust && chi2 > chi2Threshold) {
996 : // std::cout << "robust fitter called... " << std::endl;
997 0 : fitterG->EvalRobust(robustFraction);
998 0 : fitterG->GetParameters(fitParam);
999 0 : }
1000 : } else {
1001 : // set parameteres to 0
1002 : Int_t nParameters = 3;
1003 0 : if (fitType == 1)
1004 : nParameters = 6;
1005 :
1006 0 : for(Int_t i = 0; i < nParameters; i++)
1007 0 : fitParam[i] = 0;
1008 : }
1009 :
1010 0 : delete fitterG;
1011 0 : }
1012 :
1013 : AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector, EPadType padType, AliTPCCalROC *oldTPCCalROC ){
1014 : /// Create ROC with global fit parameters
1015 : /// The origin of the fit function is the center of the ROC!
1016 : /// loop over all channels, write fit values into new ROC and return it
1017 : /// In 2015 different gain in medium and long pads
1018 : /// New CalROC is created for medium pads.
1019 : /// For long pads oldTPCCalROC already contains values for medium pads.
1020 :
1021 0 : if(sector<36) padType=kAll;
1022 :
1023 : Float_t dlx, dly;
1024 0 : Float_t centerPad[3] = {0};
1025 0 : Float_t localXY[3] = {0};
1026 : AliTPCCalROC * xROCfitted = NULL;
1027 0 : if(oldTPCCalROC!=0 && padType!=kAll) xROCfitted = oldTPCCalROC;
1028 0 : else xROCfitted = new AliTPCCalROC(sector);
1029 0 : AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
1030 0 : if(padType==kAll)
1031 0 : tpcROCinstance->GetPositionLocal(sector, xROCfitted->GetNrows()/2, xROCfitted->GetNPads(xROCfitted->GetNrows()/2)/2, centerPad); // calculate center of ROC
1032 0 : else if(padType==kOROCmedium)
1033 0 : tpcROCinstance->GetPositionLocal(sector, 64/2, xROCfitted->GetNPads(64/2)/2, centerPad); // calculate center of ROC medium pads
1034 0 : else if(padType==kOROClong)
1035 0 : tpcROCinstance->GetPositionLocal(sector, (xROCfitted->GetNrows()-64)/2+64, xROCfitted->GetNPads((xROCfitted->GetNrows()-64)/2+64)/2, centerPad); // calculate center of ROC long pads
1036 :
1037 : UInt_t irowMin = 0;
1038 0 : UInt_t irowMax = xROCfitted->GetNrows();
1039 0 : if(padType==kOROCmedium) irowMax = 64;
1040 0 : else if(padType==kOROClong) irowMin = 64;
1041 :
1042 : Int_t fitType = 1;
1043 0 : if (fitParam.GetNoElements() == 6) fitType = 1;
1044 : else fitType = 0;
1045 : Double_t value = 0;
1046 0 : if (fitType == 1) { // parabolic fit
1047 0 : for (UInt_t irow = irowMin; irow < irowMax; irow++) {
1048 0 : for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
1049 0 : tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
1050 0 : dlx = localXY[0] - centerPad[0];
1051 0 : dly = localXY[1] - centerPad[1];
1052 0 : value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
1053 0 : xROCfitted->SetValue(irow, ipad, value);
1054 : }
1055 : }
1056 0 : }
1057 : else { // linear fit
1058 0 : for (UInt_t irow = irowMin; irow < irowMax; irow++) {
1059 0 : for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
1060 0 : tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
1061 0 : dlx = localXY[0] - centerPad[0];
1062 0 : dly = localXY[1] - centerPad[1];
1063 0 : value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
1064 0 : xROCfitted->SetValue(irow, ipad, value);
1065 : }
1066 : }
1067 : }
1068 0 : return xROCfitted;
1069 0 : }
1070 :
|
