]> 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
Fix for the mirrored C-side TPC sectors (Marek). Removal of BYMirror option in reco...
[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");
eea36fac 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());
db3576a7 225 if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
ebe95b38 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 //
66a2e32f 594 if (fSector%36>17){
595 c.SetY(-c.GetY());
adefcaa6 596 }
508541c7 597
1c53abe2 598 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
1c53abe2 599 c.SetType(-(c.GetType()+3)); //edge clusters
600 }
601 if (fLoop==2) c.SetType(100);
eea36fac 602 if (!AcceptCluster(&c)) return;
1c53abe2 603
604 TClonesArray * arr = fRowCl->GetArray();
03fe1804 605 AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
d101caf3 606 // if (fRecoParam->DumpSignal() &&matrix ) {
607// Int_t nbins=0;
608// Float_t *graph =0;
609// if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
610// nbins = fMaxTime;
611// graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
612// }
613// AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
614// cl->SetInfo(info);
615// }
b127a65f 616 if (!fRecoParam->DumpSignal()) {
617 cl->SetInfo(0);
618 }
eea36fac 619
620 if (AliTPCReconstructor::StreamLevel()>1) {
621 (*fDebugStreamer)<<"Clusters"<<
622 "Cl.="<<cl<<
623 "\n";
624 }
1c53abe2 625
626 fNcluster++;
627}
628
629
630//_____________________________________________________________________________
f8aae377 631void AliTPCclustererMI::Digits2Clusters()
1c53abe2 632{
633 //-----------------------------------------------------------------
634 // This is a simple cluster finder.
635 //-----------------------------------------------------------------
1c53abe2 636
f8aae377 637 if (!fInput) {
638 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 639 return;
640 }
641
13116aec 642 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 643 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1c53abe2 644 AliSimDigits digarr, *dummy=&digarr;
645 fRowDig = dummy;
646 fInput->GetBranch("Segment")->SetAddress(&dummy);
647 Stat_t nentries = fInput->GetEntries();
648
daac2a58 649 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 650
1c53abe2 651 Int_t nclusters = 0;
13116aec 652
1c53abe2 653 for (Int_t n=0; n<nentries; n++) {
654 fInput->GetEvent(n);
508541c7 655 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 656 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
657 continue;
658 }
508541c7 659 Int_t row = fRow;
13116aec 660 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 661 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 662 //
1c53abe2 663
ebe95b38 664 fRowCl->SetID(digarr.GetID());
665 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
f8aae377 666 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 667
668
f8aae377 669 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 670 fZWidth = fParam->GetZWidth();
671 if (fSector < kNIS) {
f8aae377 672 fMaxPad = fParam->GetNPadsLow(row);
1c53abe2 673 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 674 fPadLength = fParam->GetPadPitchLength(fSector,row);
675 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 676 } else {
f8aae377 677 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 678 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 679 fPadLength = fParam->GetPadPitchLength(fSector,row);
680 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 681 }
682
683
684 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 685 fBins =new Float_t[fMaxBin];
5b33a7f5 686 fSigBins =new Int_t[fMaxBin];
687 fNSigBins = 0;
13116aec 688 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 689
690 if (digarr.First()) //MI change
691 do {
13116aec 692 Float_t dig=digarr.CurrentDigit();
f8aae377 693 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 694 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 695 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 696 Int_t bin = i*fMaxTime+j;
697 fBins[bin]=dig/gain;
698 fSigBins[fNSigBins++]=bin;
1c53abe2 699 } while (digarr.Next());
700 digarr.ExpandTrackBuffer();
701
940ed1f0 702 FindClusters(noiseROC);
ebe95b38 703 FillRow();
d101caf3 704 fRowCl->GetArray()->Clear();
88cb7938 705 nclusters+=fNcluster;
5b33a7f5 706 delete[] fBins;
707 delete[] fSigBins;
88cb7938 708 }
f8aae377 709
19dd5b2f 710 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
f8aae377 711}
712
713void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
714{
715//-----------------------------------------------------------------
22c352f8 716// This is a cluster finder for the TPC raw data.
717// The method assumes NO ordering of the altro channels.
718// The pedestal subtraction can be switched on and off
719// using an option of the TPC reconstructor
f8aae377 720//-----------------------------------------------------------------
721
f8aae377 722
f8aae377 723 fRowDig = NULL;
adefcaa6 724 AliTPCROC * roc = AliTPCROC::Instance();
22c352f8 725 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 726 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
727 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
d6834f5f 728 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
729 //
730 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 731 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
732 if (fEventHeader){
733 fTimeStamp = fEventHeader->Get("Timestamp");
734 fEventType = fEventHeader->Get("Type");
735 }
736
22c352f8 737
f8aae377 738 Int_t nclusters = 0;
88cb7938 739
daac2a58 740 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 741 const Int_t kNIS = fParam->GetNInnerSector();
742 const Int_t kNOS = fParam->GetNOuterSector();
743 const Int_t kNS = kNIS + kNOS;
744 fZWidth = fParam->GetZWidth();
745 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 746 //
747 //alocate memory for sector - maximal case
748 //
22c352f8 749 Float_t** allBins = NULL;
5b33a7f5 750 Int_t** allSigBins = NULL;
751 Int_t* allNSigBins = NULL;
adefcaa6 752 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
753 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
754 allBins = new Float_t*[nRowsMax];
5b33a7f5 755 allSigBins = new Int_t*[nRowsMax];
756 allNSigBins = new Int_t[nRowsMax];
adefcaa6 757 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
758 //
759 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
760 allBins[iRow] = new Float_t[maxBin];
adefcaa6 761 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 762 allSigBins[iRow] = new Int_t[maxBin];
763 allNSigBins[iRow]=0;
adefcaa6 764 }
765 //
22c352f8 766 // Loop over sectors
adefcaa6 767 //
22c352f8 768 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 769
940ed1f0 770 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
771 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
772 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
b24bd339 773 //check the presence of the calibration
774 if (!noiseROC ||!pedestalROC ) {
775 AliError(Form("Missing calibration per sector\t%d\n",fSector));
776 continue;
777 }
22c352f8 778 Int_t nRows = 0;
779 Int_t nDDLs = 0, indexDDL = 0;
780 if (fSector < kNIS) {
781 nRows = fParam->GetNRowLow();
782 fSign = (fSector < kNIS/2) ? 1 : -1;
783 nDDLs = 2;
784 indexDDL = fSector * 2;
785 }
786 else {
787 nRows = fParam->GetNRowUp();
788 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
789 nDDLs = 4;
790 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
791 }
792
22c352f8 793 for (Int_t iRow = 0; iRow < nRows; iRow++) {
794 Int_t maxPad;
795 if (fSector < kNIS)
796 maxPad = fParam->GetNPadsLow(iRow);
797 else
798 maxPad = fParam->GetNPadsUp(iRow);
799
800 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
22c352f8 801 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 802 allNSigBins[iRow] = 0;
22c352f8 803 }
804
805 // Loas the raw data for corresponding DDLs
806 rawReader->Reset();
362c9d61 807 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
adefcaa6 808 Int_t digCounter=0;
22c352f8 809 // Begin loop over altro data
adefcaa6 810 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
811 Float_t gain =1;
812 Int_t lastPad=-1;
22c352f8 813 while (input.Next()) {
22c352f8 814 if (input.GetSector() != fSector)
815 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
816
22c352f8 817
818 Int_t iRow = input.GetRow();
9546a7ff 819 if (iRow < 0){
820 continue;
821 }
822
3f0af4ba 823 if (iRow < 0 || iRow >= nRows){
824 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 825 iRow, 0, nRows -1));
3f0af4ba 826 continue;
827 }
adefcaa6 828 //pad
829 Int_t iPad = input.GetPad();
3f0af4ba 830 if (iPad < 0 || iPad >= nPadsMax) {
831 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 832 iPad, 0, nPadsMax-1));
3f0af4ba 833 continue;
834 }
adefcaa6 835 if (iPad!=lastPad){
836 gain = gainROC->GetValue(iRow,iPad);
837 lastPad = iPad;
838 }
839 iPad+=3;
840 //time
841 Int_t iTimeBin = input.GetTime();
daac2a58 842 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
843 continue;
22c352f8 844 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 845 iTimeBin, 0, iTimeBin -1));
daac2a58 846 }
adefcaa6 847 iTimeBin+=3;
3f0af4ba 848
adefcaa6 849 //signal
22c352f8 850 Float_t signal = input.GetSignal();
adefcaa6 851 if (!calcPedestal && signal <= zeroSup) continue;
9546a7ff 852
940ed1f0 853 if (!calcPedestal) {
5b33a7f5 854 Int_t bin = iPad*fMaxTime+iTimeBin;
855 allBins[iRow][bin] = signal/gain;
856 allSigBins[iRow][allNSigBins[iRow]++] = bin;
940ed1f0 857 }else{
858 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
859 }
aba7fdfc 860 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 861
862 // Temporary
863 digCounter++;
22c352f8 864 } // End of the loop over altro data
adefcaa6 865 //
866 //
aba7fdfc 867 //
868 //
22c352f8 869 // Now loop over rows and perform pedestal subtraction
adefcaa6 870 if (digCounter==0) continue;
aba7fdfc 871 // if (calcPedestal) {
872 if (kTRUE) {
22c352f8 873 for (Int_t iRow = 0; iRow < nRows; iRow++) {
874 Int_t maxPad;
875 if (fSector < kNIS)
876 maxPad = fParam->GetNPadsLow(iRow);
877 else
878 maxPad = fParam->GetNPadsUp(iRow);
879
940ed1f0 880 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
aba7fdfc 881 //
882 // Temporary fix for data production - !!!! MARIAN
883 // The noise calibration should take mean and RMS - currently the Gaussian fit used
884 // In case of double peak - the pad should be rejected
885 //
886 // Line mean - if more than given digits over threshold - make a noise calculation
887 // and pedestal substration
888 if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
889 //
940ed1f0 890 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 891 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 892 //Float_t pedestal = TMath::Median(fMaxTime, p);
893 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 894 // calib values
895 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
896 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
897 Double_t rmsEvent = rmsCalib;
898 Double_t pedestalEvent = pedestalCalib;
899 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
900 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
901 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
902
903 //
22c352f8 904 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
5b33a7f5 905 Int_t bin = iPad*fMaxTime+iTimeBin;
906 allBins[iRow][bin] -= pedestalEvent;
7fef31a6 907 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
5b33a7f5 908 allBins[iRow][bin] = 0;
7fef31a6 909 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
5b33a7f5 910 allBins[iRow][bin] = 0;
22c352f8 911 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
5b33a7f5 912 allBins[iRow][bin] = 0;
913 if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
914 allBins[iRow][bin] = 0;
915 if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
22c352f8 916 }
917 }
f8aae377 918 }
22c352f8 919 }
eea36fac 920
921 if (AliTPCReconstructor::StreamLevel()>3) {
922 for (Int_t iRow = 0; iRow < nRows; iRow++) {
923 Int_t maxPad;
924 if (fSector < kNIS)
925 maxPad = fParam->GetNPadsLow(iRow);
926 else
927 maxPad = fParam->GetNPadsUp(iRow);
928
929 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
930 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
931 Int_t bin = iPad*fMaxTime+iTimeBin;
932 Float_t signal = allBins[iRow][bin];
933 if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
934 Double_t x[]={iRow,iPad-3,iTimeBin-3};
935 Int_t i[]={fSector};
936 AliTPCTransform trafo;
937 trafo.Transform(x,i,0,1);
938 Double_t gx[3]={x[0],x[1],x[2]};
939 trafo.RotatedGlobal2Global(fSector,gx);
940
941 (*fDebugStreamer)<<"Digits"<<
942 "sec="<<fSector<<
943 "row="<<iRow<<
944 "pad="<<iPad<<
945 "time="<<iTimeBin<<
946 "sig="<<signal<<
947 "x="<<x[0]<<
948 "y="<<x[1]<<
949 "z="<<x[2]<<
950 "gx="<<gx[0]<<
951 "gy="<<gx[1]<<
952 "gz="<<gx[2]<<
953 "\n";
954 }
955 }
956 }
957 }
958 }
959
22c352f8 960 // Now loop over rows and find clusters
961 for (fRow = 0; fRow < nRows; fRow++) {
22c352f8 962 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
ebe95b38 963 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
22c352f8 964
965 fRx = fParam->GetPadRowRadii(fSector, fRow);
966 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 967 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 968 if (fSector < kNIS)
969 fMaxPad = fParam->GetNPadsLow(fRow);
970 else
971 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 972 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 973
974 fBins = allBins[fRow];
5b33a7f5 975 fSigBins = allSigBins[fRow];
976 fNSigBins = allNSigBins[fRow];
22c352f8 977
940ed1f0 978 FindClusters(noiseROC);
ebe95b38 979 FillRow();
d101caf3 980 fRowCl->GetArray()->Clear();
22c352f8 981 nclusters += fNcluster;
982 } // End of loop to find clusters
22c352f8 983 } // End of loop over sectors
adefcaa6 984
985 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
986 delete [] allBins[iRow];
5b33a7f5 987 delete [] allSigBins[iRow];
adefcaa6 988 }
989 delete [] allBins;
5b33a7f5 990 delete [] allSigBins;
991 delete [] allNSigBins;
adefcaa6 992
9546a7ff 993 if (rawReader->GetEventId() && fOutput ){
994 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
9546a7ff 995 }
ebe95b38 996
db3576a7 997 if(rawReader->GetEventId()) {
998 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), nclusters);
999 }
f8aae377 1000}
1001
940ed1f0 1002void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 1003{
940ed1f0 1004
1005 //
1006 // add virtual charge at the edge
1007 //
7041b196 1008 Double_t kMaxDumpSize = 500000;
ebe95b38 1009 if (!fOutput) {
1010 fBDumpSignal =kFALSE;
1011 }else{
1012 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1013 }
1014
f8aae377 1015 fNcluster=0;
f8aae377 1016 fLoop=1;
a1ec4d07 1017 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 1018 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1019 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1020 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1021 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1022 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1023 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 1024 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1025 Int_t i = fSigBins[iSig];
1026 if (i%fMaxTime<=crtime) continue;
1027 Float_t *b = &fBins[i];
940ed1f0 1028 //absolute custs
03fe1804 1029 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1030 //
1031 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
aba7fdfc 1032 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1033 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 1034 //
1035 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1036 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1037 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1038 //
1039 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1040 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1041 // sigma cuts
1042 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1043 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1044 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1045
d101caf3 1046 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1047 Int_t dummy=0;
1048 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1049
f8aae377 1050 //}
1051 }
8569a2b0 1052}
65c045f0 1053
eea36fac 1054Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1055 //
679026e6 1056 // Currently hack to filter digital noise (15.06.2008)
eea36fac 1057 // To be parameterized in the AliTPCrecoParam
1058 // More inteligent way to be used in future
679026e6 1059 // Acces to the proper pedestal file needed
eea36fac 1060 //
1061 if (cl->GetMax()<400) return kTRUE;
1062 Double_t ratio = cl->GetQ()/cl->GetMax();
679026e6 1063 if (cl->GetMax()>700){
1064 if ((ratio - int(ratio)>0.8)) return kFALSE;
1065 }
eea36fac 1066 if ((ratio - int(ratio)<0.95)) return kTRUE;
1067 return kFALSE;
eea36fac 1068}
1069
65c045f0 1070
940ed1f0 1071Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1072 //
1073 // process signal on given pad - + streaming of additional information in special mode
1074 //
1075 // id[0] - sector
1076 // id[1] - row
adefcaa6 1077 // id[2] - pad
1078
1079 //
1080 // ESTIMATE pedestal and the noise
1081 //
adefcaa6 1082 const Int_t kPedMax = 100;
1083 Float_t max = 0;
1084 Float_t maxPos = 0;
1085 Int_t median = -1;
1086 Int_t count0 = 0;
1087 Int_t count1 = 0;
940ed1f0 1088 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1089 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1090 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 1091 //
adefcaa6 1092 UShort_t histo[kPedMax];
f174362f 1093 //memset(histo,0,kPedMax*sizeof(UShort_t));
b842afde 1094 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
adefcaa6 1095 for (Int_t i=0; i<fMaxTime; i++){
1096 if (signal[i]<=0) continue;
940ed1f0 1097 if (signal[i]>max && i>firstBin) {
65c045f0 1098 max = signal[i];
1099 maxPos = i;
adefcaa6 1100 }
1101 if (signal[i]>kPedMax-1) continue;
1102 histo[int(signal[i]+0.5)]++;
1103 count0++;
65c045f0 1104 }
7fef31a6 1105 //
adefcaa6 1106 for (Int_t i=1; i<kPedMax; i++){
1107 if (count1<count0*0.5) median=i;
1108 count1+=histo[i];
1109 }
1110 // truncated mean
65c045f0 1111 //
7041b196 1112 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1113 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1114 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1115 //
1116 for (Int_t idelta=1; idelta<10; idelta++){
1117 if (median-idelta<=0) continue;
1118 if (median+idelta>kPedMax) continue;
1119 if (count06<0.6*count1){
1120 count06+=histo[median-idelta];
1121 mean06 +=histo[median-idelta]*(median-idelta);
1122 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1123 count06+=histo[median+idelta];
1124 mean06 +=histo[median+idelta]*(median+idelta);
1125 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1126 }
1127 if (count09<0.9*count1){
1128 count09+=histo[median-idelta];
1129 mean09 +=histo[median-idelta]*(median-idelta);
1130 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1131 count09+=histo[median+idelta];
1132 mean09 +=histo[median+idelta]*(median+idelta);
1133 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1134 }
adefcaa6 1135 if (count10<0.95*count1){
1136 count10+=histo[median-idelta];
1137 mean +=histo[median-idelta]*(median-idelta);
1138 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1139 count10+=histo[median+idelta];
1140 mean +=histo[median+idelta]*(median+idelta);
1141 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1142 }
1143 }
3f0af4ba 1144 if (count10) {
1145 mean /=count10;
1146 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1147 }
1148 if (count06) {
1149 mean06/=count06;
1150 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1151 }
1152 if (count09) {
1153 mean09/=count09;
1154 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1155 }
940ed1f0 1156 rmsEvent = rms09;
adefcaa6 1157 //
940ed1f0 1158 pedestalEvent = median;
adefcaa6 1159 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1160 //
1161 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1162 //
1163 // Dump mean signal info
1164 //
1165 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1166 "TimeStamp="<<fTimeStamp<<
1167 "EventType="<<fEventType<<
adefcaa6 1168 "Sector="<<uid[0]<<
1169 "Row="<<uid[1]<<
1170 "Pad="<<uid[2]<<
1171 "Max="<<max<<
1172 "MaxPos="<<maxPos<<
1173 //
1174 "Median="<<median<<
1175 "Mean="<<mean<<
1176 "RMS="<<rms<<
1177 "Mean06="<<mean06<<
1178 "RMS06="<<rms06<<
1179 "Mean09="<<mean09<<
1180 "RMS09="<<rms09<<
940ed1f0 1181 "RMSCalib="<<rmsCalib<<
1182 "PedCalib="<<pedestalCalib<<
adefcaa6 1183 "\n";
1184 //
1185 // fill pedestal histogram
1186 //
65c045f0 1187 //
adefcaa6 1188 //
1189 //
1190 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1191 Float_t *dsignal = new Float_t[nchannels];
1192 Float_t *dtime = new Float_t[nchannels];
1193 for (Int_t i=0; i<nchannels; i++){
1194 dtime[i] = i;
1195 dsignal[i] = signal[i];
1196 }
03fe1804 1197
11441748 1198 TGraph * graph=0;
7fef31a6 1199 //
a525cc34 1200 // Big signals dumping
1201 //
1202 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1203 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1204 "TimeStamp="<<fTimeStamp<<
1205 "EventType="<<fEventType<<
1206 "Sector="<<uid[0]<<
1207 "Row="<<uid[1]<<
1208 "Pad="<<uid[2]<<
1209 "Graph="<<graph<<
1210 "Max="<<max<<
1211 "MaxPos="<<maxPos<<
1212 //
1213 "Median="<<median<<
1214 "Mean="<<mean<<
1215 "RMS="<<rms<<
1216 "Mean06="<<mean06<<
1217 "RMS06="<<rms06<<
1218 "Mean09="<<mean09<<
1219 "RMS09="<<rms09<<
1220 "\n";
1221 delete graph;
7fef31a6 1222
65c045f0 1223 delete [] dsignal;
1224 delete [] dtime;
940ed1f0 1225 if (rms06>fRecoParam->GetMaxNoise()) {
1226 pedestalEvent+=1024.;
1227 return 1024+median; // sign noisy channel in debug mode
1228 }
65c045f0 1229 return median;
1230}
1231
1232
1233
03fe1804 1234
03fe1804 1235