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