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