]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCclustererMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
No optimization with gcc 4.3.0
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
CommitLineData
1c53abe2 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
88cb7938 16/* $Id$ */
1c53abe2 17
18//-------------------------------------------------------
19// Implementation of the TPC clusterer
20//
ebe95b38 21// 1. The Input data for reconstruction - Options
22// 1.a Simulated data - TTree - invoked Digits2Clusters()
23// 1.b Raw data - Digits2Clusters(AliRawReader* rawReader);
24//
25// 2. The Output data
26// 2.a TTree with clusters - if SetOutput(TTree * tree) invoked
27// 2.b TObjArray - Faster option for HLT
28//
29// 3. Reconstruction setup
30// see AliTPCRecoParam for list of parameters
31// The reconstruction parameterization taken from the
32// AliTPCReconstructor::GetRecoParam()
33// Possible to setup it in reconstruction macro AliTPCReconstructor::SetRecoParam(...)
34//
35//
36//
1c53abe2 37// Origin: Marian Ivanov
38//-------------------------------------------------------
39
a1e17193 40#include "Riostream.h"
41#include <TF1.h>
1c53abe2 42#include <TFile.h>
a1e17193 43#include <TGraph.h>
44#include <TH1F.h>
45#include <TObjArray.h>
46#include <TRandom.h>
47#include <TTree.h>
48#include <TTreeStream.h>
65c045f0 49
1c53abe2 50#include "AliDigits.h"
a1e17193 51#include "AliLoader.h"
52#include "AliLog.h"
53#include "AliMathBase.h"
ef4ad662 54#include "AliRawEventHeaderBase.h"
a1e17193 55#include "AliRawReader.h"
f8aae377 56#include "AliRunLoader.h"
a1e17193 57#include "AliSimDigits.h"
13116aec 58#include "AliTPCCalPad.h"
59#include "AliTPCCalROC.h"
a1e17193 60#include "AliTPCClustersArray.h"
61#include "AliTPCClustersRow.h"
62#include "AliTPCParam.h"
63#include "AliTPCRawStream.h"
64#include "AliTPCRecoParam.h"
65#include "AliTPCReconstructor.h"
66#include "AliTPCcalibDB.h"
67#include "AliTPCclusterInfo.h"
68#include "AliTPCclusterMI.h"
022ee144 69#include "AliTPCTransform.h"
a1e17193 70#include "AliTPCclustererMI.h"
13116aec 71
1c53abe2 72ClassImp(AliTPCclustererMI)
73
74
75
2f32844a 76AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
77 fBins(0),
5b33a7f5 78 fSigBins(0),
79 fNSigBins(0),
2f32844a 80 fLoop(0),
81 fMaxBin(0),
82 fMaxTime(0),
83 fMaxPad(0),
84 fSector(-1),
85 fRow(-1),
86 fSign(0),
87 fRx(0),
88 fPadWidth(0),
89 fPadLength(0),
90 fZWidth(0),
91 fPedSubtraction(kFALSE),
ef4ad662 92 fEventHeader(0),
93 fTimeStamp(0),
94 fEventType(0),
2f32844a 95 fInput(0),
96 fOutput(0),
ebe95b38 97 fOutputArray(0),
2f32844a 98 fRowCl(0),
99 fRowDig(0),
100 fParam(0),
101 fNcluster(0),
2f32844a 102 fDebugStreamer(0),
03fe1804 103 fRecoParam(0),
a525cc34 104 fBDumpSignal(kFALSE)
1c53abe2 105{
97f127bb 106 //
107 // COSNTRUCTOR
108 // param - tpc parameters for given file
109 // recoparam - reconstruction parameters
110 //
1c53abe2 111 fInput =0;
f8aae377 112 fParam = par;
194b0609 113 if (recoParam) {
114 fRecoParam = recoParam;
115 }else{
116 //set default parameters if not specified
117 fRecoParam = AliTPCReconstructor::GetRecoParam();
118 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
119 }
65c045f0 120 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
03fe1804 121 Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
d101caf3 122 fRowCl= new AliTPCClustersRow();
123 fRowCl->SetClass("AliTPCclusterMI");
124 fRowCl->SetArray(1);
125
65c045f0 126}
e046d791 127//______________________________________________________________
128AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
129 :TObject(param),
130 fBins(0),
5b33a7f5 131 fSigBins(0),
132 fNSigBins(0),
e046d791 133 fLoop(0),
134 fMaxBin(0),
135 fMaxTime(0),
136 fMaxPad(0),
137 fSector(-1),
138 fRow(-1),
139 fSign(0),
140 fRx(0),
141 fPadWidth(0),
142 fPadLength(0),
143 fZWidth(0),
144 fPedSubtraction(kFALSE),
e046d791 145 fEventHeader(0),
146 fTimeStamp(0),
147 fEventType(0),
148 fInput(0),
149 fOutput(0),
ebe95b38 150 fOutputArray(0),
e046d791 151 fRowCl(0),
152 fRowDig(0),
153 fParam(0),
154 fNcluster(0),
e046d791 155 fDebugStreamer(0),
5b33a7f5 156 fRecoParam(0),
a525cc34 157 fBDumpSignal(kFALSE)
e046d791 158{
159 //
160 // dummy
161 //
162 fMaxBin = param.fMaxBin;
163}
164//______________________________________________________________
165AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
166{
167 //
168 // assignment operator - dummy
169 //
170 fMaxBin=param.fMaxBin;
171 return (*this);
172}
173//______________________________________________________________
65c045f0 174AliTPCclustererMI::~AliTPCclustererMI(){
ebe95b38 175 //
176 //
177 //
65c045f0 178 if (fDebugStreamer) delete fDebugStreamer;
ebe95b38 179 if (fOutputArray){
af86c1fd 180 //fOutputArray->Delete();
ebe95b38 181 delete fOutputArray;
182 }
1c53abe2 183}
22c352f8 184
1c53abe2 185void AliTPCclustererMI::SetInput(TTree * tree)
186{
187 //
188 // set input tree with digits
189 //
190 fInput = tree;
191 if (!fInput->GetBranch("Segment")){
192 cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
193 fInput=0;
194 return;
195 }
196}
197
198void AliTPCclustererMI::SetOutput(TTree * tree)
199{
200 //
ebe95b38 201 // Set the output tree
202 // If not set the ObjArray used - Option for HLT
1c53abe2 203 //
ebe95b38 204 if (!tree) return;
205 fOutput= tree;
1c53abe2 206 AliTPCClustersRow clrow;
207 AliTPCClustersRow *pclrow=&clrow;
208 clrow.SetClass("AliTPCclusterMI");
209 clrow.SetArray(1); // to make Clones array
210 fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
211}
212
213
ebe95b38 214void AliTPCclustererMI::FillRow(){
215 //
216 // fill the output container -
217 // 2 Options possible
218 // Tree
219 // TObjArray
220 //
221 if (fOutput) fOutput->Fill();
222 if (!fOutput){
223 //
af86c1fd 224 if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
ebe95b38 225 if (fRowCl) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
226 }
227}
228
1c53abe2 229Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
230 // sigma y2 = in digits - we don't know the angle
753797ce 231 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 232 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
233 (fPadWidth*fPadWidth);
234 Float_t sres = 0.25;
235 Float_t res = sd2+sres;
236 return res;
237}
238
239
240Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){
241 //sigma z2 = in digits - angle estimated supposing vertex constraint
753797ce 242 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 243 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
a1ec4d07 244 Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
1c53abe2 245 angular*=angular;
246 angular/=12.;
247 Float_t sres = fParam->GetZSigma()/fZWidth;
248 sres *=sres;
249 Float_t res = angular +sd2+sres;
250 return res;
251}
252
13116aec 253void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
1c53abe2 254AliTPCclusterMI &c)
255{
97f127bb 256 //
257 // k - Make cluster at position k
258 // bins - 2 D array of signals mapped to 1 dimensional array -
259 // max - the number of time bins er one dimension
260 // c - refernce to cluster to be filled
261 //
1c53abe2 262 Int_t i0=k/max; //central pad
263 Int_t j0=k%max; //central time bin
264
265 // set pointers to data
266 //Int_t dummy[5] ={0,0,0,0,0};
13116aec 267 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
1c53abe2 268 for (Int_t di=-2;di<=2;di++){
269 matrix[di+2] = &bins[k+di*max];
1c53abe2 270 }
271 //build matrix with virtual charge
272 Float_t sigmay2= GetSigmaY2(j0);
273 Float_t sigmaz2= GetSigmaZ2(j0);
274
275 Float_t vmatrix[5][5];
276 vmatrix[2][2] = matrix[2][0];
277 c.SetType(0);
6c024a0e 278 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
1c53abe2 279 for (Int_t di =-1;di <=1;di++)
280 for (Int_t dj =-1;dj <=1;dj++){
281 Float_t amp = matrix[di+2][dj];
282 if ( (amp<2) && (fLoop<2)){
283 // if under threshold - calculate virtual charge
284 Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
285 amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
286 if (amp>2) amp = 2;
287 vmatrix[2+di][2+dj]=amp;
288 vmatrix[2+2*di][2+2*dj]=0;
289 if ( (di*dj)!=0){
290 //DIAGONAL ELEMENTS
291 vmatrix[2+2*di][2+dj] =0;
292 vmatrix[2+di][2+2*dj] =0;
293 }
294 continue;
295 }
296 if (amp<4){
297 //if small amplitude - below 2 x threshold - don't consider other one
298 vmatrix[2+di][2+dj]=amp;
299 vmatrix[2+2*di][2+2*dj]=0; // don't take to the account next bin
300 if ( (di*dj)!=0){
301 //DIAGONAL ELEMENTS
302 vmatrix[2+2*di][2+dj] =0;
303 vmatrix[2+di][2+2*dj] =0;
304 }
305 continue;
306 }
307 //if bigger then take everything
308 vmatrix[2+di][2+dj]=amp;
309 vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
310 if ( (di*dj)!=0){
311 //DIAGONAL ELEMENTS
312 vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
313 vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
314 }
315 }
316
317
318
319 Float_t sumw=0;
320 Float_t sumiw=0;
321 Float_t sumi2w=0;
322 Float_t sumjw=0;
323 Float_t sumj2w=0;
324 //
325 for (Int_t i=-2;i<=2;i++)
326 for (Int_t j=-2;j<=2;j++){
327 Float_t amp = vmatrix[i+2][j+2];
328
329 sumw += amp;
330 sumiw += i*amp;
331 sumi2w += i*i*amp;
332 sumjw += j*amp;
333 sumj2w += j*j*amp;
334 }
335 //
336 Float_t meani = sumiw/sumw;
337 Float_t mi2 = sumi2w/sumw-meani*meani;
338 Float_t meanj = sumjw/sumw;
339 Float_t mj2 = sumj2w/sumw-meanj*meanj;
340 //
341 Float_t ry = mi2/sigmay2;
342 Float_t rz = mj2/sigmaz2;
343
344 //
345 if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
194b0609 346 if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) {
347 //
348 //if cluster looks like expected or Unfolding not switched on
349 //standard COG is used
1c53abe2 350 //+1.2 deviation from expected sigma accepted
351 // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
352
353 meani +=i0;
354 meanj +=j0;
355 //set cluster parameters
356 c.SetQ(sumw);
022ee144 357 c.SetPad(meani-2.5);
5bd7ed29 358 c.SetTimeBin(meanj-3);
1c53abe2 359 c.SetSigmaY2(mi2);
360 c.SetSigmaZ2(mj2);
022ee144 361 c.SetType(0);
03fe1804 362 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 363 return;
364 }
365 //
366 //unfolding when neccessary
367 //
368
13116aec 369 Float_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3
370 Float_t dummy[7]={0,0,0,0,0,0};
1c53abe2 371 for (Int_t di=-3;di<=3;di++){
372 matrix2[di+3] = &bins[k+di*max];
373 if ((k+di*max)<3) matrix2[di+3] = &dummy[3];
374 if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3];
375 }
376 Float_t vmatrix2[5][5];
377 Float_t sumu;
378 Float_t overlap;
379 UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
380 //
381 // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
382 meani +=i0;
383 meanj +=j0;
384 //set cluster parameters
385 c.SetQ(sumu);
022ee144 386 c.SetPad(meani-2.5);
387 c.SetTimeBin(meanj-3);
1c53abe2 388 c.SetSigmaY2(mi2);
389 c.SetSigmaZ2(mj2);
390 c.SetType(Char_t(overlap)+1);
03fe1804 391 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 392
393 //unfolding 2
394 meani-=i0;
395 meanj-=j0;
396 if (gDebug>4)
397 printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
398}
399
400
401
13116aec 402void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
1c53abe2 403 Float_t & sumu, Float_t & overlap )
404{
405 //
406 //unfold cluster from input matrix
407 //data corresponding to cluster writen in recmatrix
408 //output meani and meanj
409
410 //take separatelly y and z
411
412 Float_t sum3i[7] = {0,0,0,0,0,0,0};
413 Float_t sum3j[7] = {0,0,0,0,0,0,0};
414
415 for (Int_t k =0;k<7;k++)
416 for (Int_t l = -1; l<=1;l++){
417 sum3i[k]+=matrix2[k][l];
418 sum3j[k]+=matrix2[l+3][k-3];
419 }
420 Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
421 //
422 //unfold y
423 Float_t sum3wi = 0; //charge minus overlap
424 Float_t sum3wio = 0; //full charge
425 Float_t sum3iw = 0; //sum for mean value
426 for (Int_t dk=-1;dk<=1;dk++){
427 sum3wio+=sum3i[dk+3];
428 if (dk==0){
429 sum3wi+=sum3i[dk+3];
430 }
431 else{
432 Float_t ratio =1;
433 if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
434 sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
435 Float_t xm2 = sum3i[-dk+3];
436 Float_t xm1 = sum3i[+3];
437 Float_t x1 = sum3i[2*dk+3];
438 Float_t x2 = sum3i[3*dk+3];
cc5e9db0 439 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
440 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 441 ratio = w11/(w11+w12);
442 for (Int_t dl=-1;dl<=1;dl++)
443 mratio[dk+1][dl+1] *= ratio;
444 }
445 Float_t amp = sum3i[dk+3]*ratio;
446 sum3wi+=amp;
447 sum3iw+= dk*amp;
448 }
449 }
450 meani = sum3iw/sum3wi;
451 Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
452
453
454
455 //unfold z
456 Float_t sum3wj = 0; //charge minus overlap
457 Float_t sum3wjo = 0; //full charge
458 Float_t sum3jw = 0; //sum for mean value
459 for (Int_t dk=-1;dk<=1;dk++){
460 sum3wjo+=sum3j[dk+3];
461 if (dk==0){
462 sum3wj+=sum3j[dk+3];
463 }
464 else{
465 Float_t ratio =1;
466 if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
467 (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
468 Float_t xm2 = sum3j[-dk+3];
469 Float_t xm1 = sum3j[+3];
470 Float_t x1 = sum3j[2*dk+3];
471 Float_t x2 = sum3j[3*dk+3];
cc5e9db0 472 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
473 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 474 ratio = w11/(w11+w12);
475 for (Int_t dl=-1;dl<=1;dl++)
476 mratio[dl+1][dk+1] *= ratio;
477 }
478 Float_t amp = sum3j[dk+3]*ratio;
479 sum3wj+=amp;
480 sum3jw+= dk*amp;
481 }
482 }
483 meanj = sum3jw/sum3wj;
484 Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
485 overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
486 sumu = (sum3wj+sum3wi)/2.;
487
488 if (overlap ==3) {
489 //if not overlap detected remove everything
490 for (Int_t di =-2; di<=2;di++)
491 for (Int_t dj =-2; dj<=2;dj++){
492 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
493 }
494 }
495 else{
496 for (Int_t di =-1; di<=1;di++)
497 for (Int_t dj =-1; dj<=1;dj++){
498 Float_t ratio =1;
499 if (mratio[di+1][dj+1]==1){
500 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
501 if (TMath::Abs(di)+TMath::Abs(dj)>1){
502 recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
503 recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
504 }
505 recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
506 }
507 else
508 {
509 //if we have overlap in direction
510 recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
511 if (TMath::Abs(di)+TMath::Abs(dj)>1){
cc5e9db0 512 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
1c53abe2 513 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
514 //
cc5e9db0 515 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
1c53abe2 516 recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
517 }
518 else{
519 ratio = recmatrix[di+2][dj+2]/matrix2[3][0];
520 recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
521 }
522 }
523 }
524 }
525 if (gDebug>4)
526 printf("%f\n", recmatrix[2][2]);
527
528}
529
530Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
531{
532 //
533 // estimate max
534 Float_t sumteor= 0;
535 Float_t sumamp = 0;
536
537 for (Int_t di = -1;di<=1;di++)
538 for (Int_t dj = -1;dj<=1;dj++){
539 if (vmatrix[2+di][2+dj]>2){
540 Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
541 sumteor += teor*vmatrix[2+di][2+dj];
542 sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
543 }
544 }
545 Float_t max = sumamp/sumteor;
546 return max;
547}
548
03fe1804 549void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
1c53abe2 550 //
1c53abe2 551 //
022ee144 552 // Transform cluster to the rotated global coordinata
553 // Assign labels to the cluster
554 // add the cluster to the array
555 // for more details - See AliTPCTranform::Transform(x,i,0,1)
556 Float_t meani = c.GetPad();
557 Float_t meanj = c.GetTimeBin();
1c53abe2 558
022ee144 559 Int_t ki = TMath::Nint(meani);
1c53abe2 560 if (ki<0) ki=0;
561 if (ki>=fMaxPad) ki = fMaxPad-1;
022ee144 562 Int_t kj = TMath::Nint(meanj);
1c53abe2 563 if (kj<0) kj=0;
564 if (kj>=fMaxTime-3) kj=fMaxTime-4;
022ee144 565 // ki and kj shifted as integers coordinata
f8aae377 566 if (fRowDig) {
567 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
568 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
569 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
570 }
1c53abe2 571
022ee144 572 c.SetRow(fRow);
573 c.SetDetector(fSector);
1c53abe2 574 Float_t s2 = c.GetSigmaY2();
575 Float_t w=fParam->GetPadPitchWidth(fSector);
1c53abe2 576 c.SetSigmaY2(s2*w*w);
577 s2 = c.GetSigmaZ2();
022ee144 578 c.SetSigmaZ2(s2*fZWidth*fZWidth);
579 //
580 //
581 //
582 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
583 if (!transform) {
584 AliFatal("Tranformations not in calibDB");
585 }
586 Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
587 Int_t i[1]={fSector};
588 transform->Transform(x,i,0,1);
589 c.SetX(x[0]);
590 c.SetY(x[1]);
591 c.SetZ(x[2]);
592 //
593 //
adefcaa6 594 if (!fRecoParam->GetBYMirror()){
595 if (fSector%36>17){
022ee144 596 c.SetY(-c.GetY());
adefcaa6 597 }
598 }
508541c7 599
1c53abe2 600 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
1c53abe2 601 c.SetType(-(c.GetType()+3)); //edge clusters
602 }
603 if (fLoop==2) c.SetType(100);
604
605 TClonesArray * arr = fRowCl->GetArray();
03fe1804 606 AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
d101caf3 607 // if (fRecoParam->DumpSignal() &&matrix ) {
608// Int_t nbins=0;
609// Float_t *graph =0;
610// if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
611// nbins = fMaxTime;
612// graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
613// }
614// AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
615// cl->SetInfo(info);
616// }
b127a65f 617 if (!fRecoParam->DumpSignal()) {
618 cl->SetInfo(0);
619 }
1c53abe2 620
621 fNcluster++;
622}
623
624
625//_____________________________________________________________________________
f8aae377 626void AliTPCclustererMI::Digits2Clusters()
1c53abe2 627{
628 //-----------------------------------------------------------------
629 // This is a simple cluster finder.
630 //-----------------------------------------------------------------
1c53abe2 631
f8aae377 632 if (!fInput) {
633 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 634 return;
635 }
636
13116aec 637 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 638 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1c53abe2 639 AliSimDigits digarr, *dummy=&digarr;
640 fRowDig = dummy;
641 fInput->GetBranch("Segment")->SetAddress(&dummy);
642 Stat_t nentries = fInput->GetEntries();
643
daac2a58 644 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 645
1c53abe2 646 Int_t nclusters = 0;
13116aec 647
1c53abe2 648 for (Int_t n=0; n<nentries; n++) {
649 fInput->GetEvent(n);
508541c7 650 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 651 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
652 continue;
653 }
508541c7 654 Int_t row = fRow;
13116aec 655 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 656 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 657 //
1c53abe2 658
ebe95b38 659 fRowCl->SetID(digarr.GetID());
660 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
f8aae377 661 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 662
663
f8aae377 664 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 665 fZWidth = fParam->GetZWidth();
666 if (fSector < kNIS) {
f8aae377 667 fMaxPad = fParam->GetNPadsLow(row);
1c53abe2 668 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 669 fPadLength = fParam->GetPadPitchLength(fSector,row);
670 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 671 } else {
f8aae377 672 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 673 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 674 fPadLength = fParam->GetPadPitchLength(fSector,row);
675 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 676 }
677
678
679 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 680 fBins =new Float_t[fMaxBin];
5b33a7f5 681 fSigBins =new Int_t[fMaxBin];
682 fNSigBins = 0;
13116aec 683 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 684
685 if (digarr.First()) //MI change
686 do {
13116aec 687 Float_t dig=digarr.CurrentDigit();
f8aae377 688 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 689 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 690 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 691 Int_t bin = i*fMaxTime+j;
692 fBins[bin]=dig/gain;
693 fSigBins[fNSigBins++]=bin;
1c53abe2 694 } while (digarr.Next());
695 digarr.ExpandTrackBuffer();
696
940ed1f0 697 FindClusters(noiseROC);
ebe95b38 698 FillRow();
d101caf3 699 fRowCl->GetArray()->Clear();
88cb7938 700 nclusters+=fNcluster;
5b33a7f5 701 delete[] fBins;
702 delete[] fSigBins;
88cb7938 703 }
f8aae377 704
19dd5b2f 705 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
f8aae377 706}
707
708void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
709{
710//-----------------------------------------------------------------
22c352f8 711// This is a cluster finder for the TPC raw data.
712// The method assumes NO ordering of the altro channels.
713// The pedestal subtraction can be switched on and off
714// using an option of the TPC reconstructor
f8aae377 715//-----------------------------------------------------------------
716
f8aae377 717
f8aae377 718 fRowDig = NULL;
adefcaa6 719 AliTPCROC * roc = AliTPCROC::Instance();
22c352f8 720 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 721 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
722 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
d6834f5f 723 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
724 //
725 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 726 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
727 if (fEventHeader){
728 fTimeStamp = fEventHeader->Get("Timestamp");
729 fEventType = fEventHeader->Get("Type");
730 }
731
22c352f8 732
f8aae377 733 Int_t nclusters = 0;
88cb7938 734
daac2a58 735 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 736 const Int_t kNIS = fParam->GetNInnerSector();
737 const Int_t kNOS = fParam->GetNOuterSector();
738 const Int_t kNS = kNIS + kNOS;
739 fZWidth = fParam->GetZWidth();
740 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 741 //
742 //alocate memory for sector - maximal case
743 //
22c352f8 744 Float_t** allBins = NULL;
5b33a7f5 745 Int_t** allSigBins = NULL;
746 Int_t* allNSigBins = NULL;
adefcaa6 747 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
748 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
749 allBins = new Float_t*[nRowsMax];
5b33a7f5 750 allSigBins = new Int_t*[nRowsMax];
751 allNSigBins = new Int_t[nRowsMax];
adefcaa6 752 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
753 //
754 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
755 allBins[iRow] = new Float_t[maxBin];
adefcaa6 756 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 757 allSigBins[iRow] = new Int_t[maxBin];
758 allNSigBins[iRow]=0;
adefcaa6 759 }
760 //
22c352f8 761 // Loop over sectors
adefcaa6 762 //
22c352f8 763 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 764
940ed1f0 765 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
766 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
767 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
b24bd339 768 //check the presence of the calibration
769 if (!noiseROC ||!pedestalROC ) {
770 AliError(Form("Missing calibration per sector\t%d\n",fSector));
771 continue;
772 }
22c352f8 773 Int_t nRows = 0;
774 Int_t nDDLs = 0, indexDDL = 0;
775 if (fSector < kNIS) {
776 nRows = fParam->GetNRowLow();
777 fSign = (fSector < kNIS/2) ? 1 : -1;
778 nDDLs = 2;
779 indexDDL = fSector * 2;
780 }
781 else {
782 nRows = fParam->GetNRowUp();
783 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
784 nDDLs = 4;
785 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
786 }
787
22c352f8 788 for (Int_t iRow = 0; iRow < nRows; iRow++) {
789 Int_t maxPad;
790 if (fSector < kNIS)
791 maxPad = fParam->GetNPadsLow(iRow);
792 else
793 maxPad = fParam->GetNPadsUp(iRow);
794
795 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
22c352f8 796 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 797 allNSigBins[iRow] = 0;
22c352f8 798 }
799
800 // Loas the raw data for corresponding DDLs
801 rawReader->Reset();
362c9d61 802 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
adefcaa6 803 Int_t digCounter=0;
22c352f8 804 // Begin loop over altro data
adefcaa6 805 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
806 Float_t gain =1;
807 Int_t lastPad=-1;
22c352f8 808 while (input.Next()) {
22c352f8 809 if (input.GetSector() != fSector)
810 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
811
22c352f8 812
813 Int_t iRow = input.GetRow();
3f0af4ba 814 if (iRow < 0 || iRow >= nRows){
815 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 816 iRow, 0, nRows -1));
3f0af4ba 817 continue;
818 }
adefcaa6 819 //pad
820 Int_t iPad = input.GetPad();
3f0af4ba 821 if (iPad < 0 || iPad >= nPadsMax) {
822 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 823 iPad, 0, nPadsMax-1));
3f0af4ba 824 continue;
825 }
adefcaa6 826 if (iPad!=lastPad){
827 gain = gainROC->GetValue(iRow,iPad);
828 lastPad = iPad;
829 }
830 iPad+=3;
831 //time
832 Int_t iTimeBin = input.GetTime();
daac2a58 833 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
834 continue;
22c352f8 835 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 836 iTimeBin, 0, iTimeBin -1));
daac2a58 837 }
adefcaa6 838 iTimeBin+=3;
3f0af4ba 839
adefcaa6 840 //signal
22c352f8 841 Float_t signal = input.GetSignal();
adefcaa6 842 if (!calcPedestal && signal <= zeroSup) continue;
940ed1f0 843 if (!calcPedestal) {
5b33a7f5 844 Int_t bin = iPad*fMaxTime+iTimeBin;
845 allBins[iRow][bin] = signal/gain;
846 allSigBins[iRow][allNSigBins[iRow]++] = bin;
940ed1f0 847 }else{
848 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
849 }
aba7fdfc 850 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 851
852 // Temporary
853 digCounter++;
22c352f8 854 } // End of the loop over altro data
adefcaa6 855 //
856 //
aba7fdfc 857 //
858 //
22c352f8 859 // Now loop over rows and perform pedestal subtraction
adefcaa6 860 if (digCounter==0) continue;
aba7fdfc 861 // if (calcPedestal) {
862 if (kTRUE) {
22c352f8 863 for (Int_t iRow = 0; iRow < nRows; iRow++) {
864 Int_t maxPad;
865 if (fSector < kNIS)
866 maxPad = fParam->GetNPadsLow(iRow);
867 else
868 maxPad = fParam->GetNPadsUp(iRow);
869
940ed1f0 870 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
aba7fdfc 871 //
872 // Temporary fix for data production - !!!! MARIAN
873 // The noise calibration should take mean and RMS - currently the Gaussian fit used
874 // In case of double peak - the pad should be rejected
875 //
876 // Line mean - if more than given digits over threshold - make a noise calculation
877 // and pedestal substration
878 if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
879 //
940ed1f0 880 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 881 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 882 //Float_t pedestal = TMath::Median(fMaxTime, p);
883 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 884 // calib values
885 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
886 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
887 Double_t rmsEvent = rmsCalib;
888 Double_t pedestalEvent = pedestalCalib;
889 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
890 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
891 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
892
893 //
22c352f8 894 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
5b33a7f5 895 Int_t bin = iPad*fMaxTime+iTimeBin;
896 allBins[iRow][bin] -= pedestalEvent;
7fef31a6 897 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
5b33a7f5 898 allBins[iRow][bin] = 0;
7fef31a6 899 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
5b33a7f5 900 allBins[iRow][bin] = 0;
22c352f8 901 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
5b33a7f5 902 allBins[iRow][bin] = 0;
903 if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
904 allBins[iRow][bin] = 0;
905 if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
22c352f8 906 }
907 }
f8aae377 908 }
22c352f8 909 }
22c352f8 910 // Now loop over rows and find clusters
911 for (fRow = 0; fRow < nRows; fRow++) {
22c352f8 912 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
ebe95b38 913 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
22c352f8 914
915 fRx = fParam->GetPadRowRadii(fSector, fRow);
916 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 917 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 918 if (fSector < kNIS)
919 fMaxPad = fParam->GetNPadsLow(fRow);
920 else
921 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 922 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 923
924 fBins = allBins[fRow];
5b33a7f5 925 fSigBins = allSigBins[fRow];
926 fNSigBins = allNSigBins[fRow];
22c352f8 927
940ed1f0 928 FindClusters(noiseROC);
ebe95b38 929 FillRow();
d101caf3 930 fRowCl->GetArray()->Clear();
22c352f8 931 nclusters += fNcluster;
932 } // End of loop to find clusters
22c352f8 933 } // End of loop over sectors
adefcaa6 934
935 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
936 delete [] allBins[iRow];
5b33a7f5 937 delete [] allSigBins[iRow];
adefcaa6 938 }
939 delete [] allBins;
5b33a7f5 940 delete [] allSigBins;
941 delete [] allNSigBins;
adefcaa6 942
d09996c4 943// if (rawReader->GetEventId() && fOutput ){
944// Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
945// }else{
946// Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), nclusters);
ebe95b38 947
d09996c4 948// }
ebe95b38 949
f8aae377 950}
951
940ed1f0 952void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 953{
940ed1f0 954
955 //
956 // add virtual charge at the edge
957 //
7041b196 958 Double_t kMaxDumpSize = 500000;
ebe95b38 959 if (!fOutput) {
960 fBDumpSignal =kFALSE;
961 }else{
962 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
963 }
964
f8aae377 965 fNcluster=0;
f8aae377 966 fLoop=1;
a1ec4d07 967 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 968 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
969 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
970 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
971 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
972 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
973 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 974 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
975 Int_t i = fSigBins[iSig];
976 if (i%fMaxTime<=crtime) continue;
977 Float_t *b = &fBins[i];
940ed1f0 978 //absolute custs
03fe1804 979 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
980 //
981 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
aba7fdfc 982 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
983 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 984 //
985 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 986 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 987 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 988 //
989 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 990 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 991 // sigma cuts
992 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
993 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
994 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
995
d101caf3 996 AliTPCclusterMI c; // default cosntruction without info
f8aae377 997 Int_t dummy=0;
998 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 999
f8aae377 1000 //}
1001 }
8569a2b0 1002}
65c045f0 1003
1004
940ed1f0 1005Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1006 //
1007 // process signal on given pad - + streaming of additional information in special mode
1008 //
1009 // id[0] - sector
1010 // id[1] - row
adefcaa6 1011 // id[2] - pad
1012
1013 //
1014 // ESTIMATE pedestal and the noise
1015 //
adefcaa6 1016 const Int_t kPedMax = 100;
7041b196 1017 Double_t kMaxDebugSize = 5000000.;
adefcaa6 1018 Float_t max = 0;
1019 Float_t maxPos = 0;
1020 Int_t median = -1;
1021 Int_t count0 = 0;
1022 Int_t count1 = 0;
940ed1f0 1023 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1024 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1025 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 1026 //
adefcaa6 1027 UShort_t histo[kPedMax];
1028 memset(histo,0,kPedMax*sizeof(UShort_t));
1029 for (Int_t i=0; i<fMaxTime; i++){
1030 if (signal[i]<=0) continue;
940ed1f0 1031 if (signal[i]>max && i>firstBin) {
65c045f0 1032 max = signal[i];
1033 maxPos = i;
adefcaa6 1034 }
1035 if (signal[i]>kPedMax-1) continue;
1036 histo[int(signal[i]+0.5)]++;
1037 count0++;
65c045f0 1038 }
7fef31a6 1039 //
adefcaa6 1040 for (Int_t i=1; i<kPedMax; i++){
1041 if (count1<count0*0.5) median=i;
1042 count1+=histo[i];
1043 }
1044 // truncated mean
65c045f0 1045 //
7041b196 1046 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1047 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1048 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1049 //
1050 for (Int_t idelta=1; idelta<10; idelta++){
1051 if (median-idelta<=0) continue;
1052 if (median+idelta>kPedMax) continue;
1053 if (count06<0.6*count1){
1054 count06+=histo[median-idelta];
1055 mean06 +=histo[median-idelta]*(median-idelta);
1056 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1057 count06+=histo[median+idelta];
1058 mean06 +=histo[median+idelta]*(median+idelta);
1059 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1060 }
1061 if (count09<0.9*count1){
1062 count09+=histo[median-idelta];
1063 mean09 +=histo[median-idelta]*(median-idelta);
1064 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1065 count09+=histo[median+idelta];
1066 mean09 +=histo[median+idelta]*(median+idelta);
1067 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1068 }
adefcaa6 1069 if (count10<0.95*count1){
1070 count10+=histo[median-idelta];
1071 mean +=histo[median-idelta]*(median-idelta);
1072 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1073 count10+=histo[median+idelta];
1074 mean +=histo[median+idelta]*(median+idelta);
1075 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1076 }
1077 }
3f0af4ba 1078 if (count10) {
1079 mean /=count10;
1080 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1081 }
1082 if (count06) {
1083 mean06/=count06;
1084 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1085 }
1086 if (count09) {
1087 mean09/=count09;
1088 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1089 }
940ed1f0 1090 rmsEvent = rms09;
adefcaa6 1091 //
940ed1f0 1092 pedestalEvent = median;
adefcaa6 1093 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1094 //
1095 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1096 //
1097 // Dump mean signal info
1098 //
1099 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1100 "TimeStamp="<<fTimeStamp<<
1101 "EventType="<<fEventType<<
adefcaa6 1102 "Sector="<<uid[0]<<
1103 "Row="<<uid[1]<<
1104 "Pad="<<uid[2]<<
1105 "Max="<<max<<
1106 "MaxPos="<<maxPos<<
1107 //
1108 "Median="<<median<<
1109 "Mean="<<mean<<
1110 "RMS="<<rms<<
1111 "Mean06="<<mean06<<
1112 "RMS06="<<rms06<<
1113 "Mean09="<<mean09<<
1114 "RMS09="<<rms09<<
940ed1f0 1115 "RMSCalib="<<rmsCalib<<
1116 "PedCalib="<<pedestalCalib<<
adefcaa6 1117 "\n";
1118 //
1119 // fill pedestal histogram
1120 //
1121 AliTPCROC * roc = AliTPCROC::Instance();
a525cc34 1122
65c045f0 1123 //
adefcaa6 1124 //
1125 //
1126 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1127 Float_t *dsignal = new Float_t[nchannels];
1128 Float_t *dtime = new Float_t[nchannels];
1129 for (Int_t i=0; i<nchannels; i++){
1130 dtime[i] = i;
1131 dsignal[i] = signal[i];
1132 }
03fe1804 1133
65c045f0 1134 TGraph * graph;
7fef31a6 1135 //
a525cc34 1136 // Big signals dumping
1137 //
1138 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1139 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1140 "TimeStamp="<<fTimeStamp<<
1141 "EventType="<<fEventType<<
1142 "Sector="<<uid[0]<<
1143 "Row="<<uid[1]<<
1144 "Pad="<<uid[2]<<
1145 "Graph="<<graph<<
1146 "Max="<<max<<
1147 "MaxPos="<<maxPos<<
1148 //
1149 "Median="<<median<<
1150 "Mean="<<mean<<
1151 "RMS="<<rms<<
1152 "Mean06="<<mean06<<
1153 "RMS06="<<rms06<<
1154 "Mean09="<<mean09<<
1155 "RMS09="<<rms09<<
1156 "\n";
1157 delete graph;
7fef31a6 1158
65c045f0 1159 delete [] dsignal;
1160 delete [] dtime;
940ed1f0 1161 if (rms06>fRecoParam->GetMaxNoise()) {
1162 pedestalEvent+=1024.;
1163 return 1024+median; // sign noisy channel in debug mode
1164 }
65c045f0 1165 return median;
1166}
1167
1168
1169
03fe1804 1170
03fe1804 1171