]> 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
syst. check that reads the MC information for the case of TPC-only
[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 }
65c045f0 125 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
eea36fac 126 // Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
d101caf3 127 fRowCl= new AliTPCClustersRow();
128 fRowCl->SetClass("AliTPCclusterMI");
129 fRowCl->SetArray(1);
130
65c045f0 131}
e046d791 132//______________________________________________________________
133AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
134 :TObject(param),
135 fBins(0),
5b33a7f5 136 fSigBins(0),
137 fNSigBins(0),
e046d791 138 fLoop(0),
139 fMaxBin(0),
140 fMaxTime(0),
141 fMaxPad(0),
142 fSector(-1),
143 fRow(-1),
144 fSign(0),
145 fRx(0),
146 fPadWidth(0),
147 fPadLength(0),
148 fZWidth(0),
149 fPedSubtraction(kFALSE),
e046d791 150 fEventHeader(0),
151 fTimeStamp(0),
152 fEventType(0),
153 fInput(0),
154 fOutput(0),
ebe95b38 155 fOutputArray(0),
98ee6d31 156 fOutputClonesArray(0),
e046d791 157 fRowCl(0),
158 fRowDig(0),
159 fParam(0),
160 fNcluster(0),
98ee6d31 161 fNclusters(0),
e046d791 162 fDebugStreamer(0),
5b33a7f5 163 fRecoParam(0),
98ee6d31 164 fBDumpSignal(kFALSE),
165 fBClonesArray(kFALSE)
e046d791 166{
167 //
168 // dummy
169 //
170 fMaxBin = param.fMaxBin;
171}
172//______________________________________________________________
173AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
174{
175 //
176 // assignment operator - dummy
177 //
178 fMaxBin=param.fMaxBin;
179 return (*this);
180}
181//______________________________________________________________
65c045f0 182AliTPCclustererMI::~AliTPCclustererMI(){
ebe95b38 183 //
184 //
185 //
65c045f0 186 if (fDebugStreamer) delete fDebugStreamer;
ebe95b38 187 if (fOutputArray){
af86c1fd 188 //fOutputArray->Delete();
ebe95b38 189 delete fOutputArray;
190 }
98ee6d31 191 if (fOutputClonesArray){
192 fOutputClonesArray->Delete();
193 delete fOutputClonesArray;
194 }
1c53abe2 195}
22c352f8 196
1c53abe2 197void AliTPCclustererMI::SetInput(TTree * tree)
198{
199 //
200 // set input tree with digits
201 //
202 fInput = tree;
203 if (!fInput->GetBranch("Segment")){
204 cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
205 fInput=0;
206 return;
207 }
208}
209
210void AliTPCclustererMI::SetOutput(TTree * tree)
211{
212 //
ebe95b38 213 // Set the output tree
214 // If not set the ObjArray used - Option for HLT
1c53abe2 215 //
ebe95b38 216 if (!tree) return;
217 fOutput= tree;
1c53abe2 218 AliTPCClustersRow clrow;
219 AliTPCClustersRow *pclrow=&clrow;
220 clrow.SetClass("AliTPCclusterMI");
221 clrow.SetArray(1); // to make Clones array
222 fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
223}
224
225
ebe95b38 226void AliTPCclustererMI::FillRow(){
227 //
228 // fill the output container -
229 // 2 Options possible
230 // Tree
231 // TObjArray
232 //
233 if (fOutput) fOutput->Fill();
98ee6d31 234 if (!fOutput && !fBClonesArray){
ebe95b38 235 //
af86c1fd 236 if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
db3576a7 237 if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
ebe95b38 238 }
239}
240
1c53abe2 241Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
242 // sigma y2 = in digits - we don't know the angle
753797ce 243 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 244 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
245 (fPadWidth*fPadWidth);
246 Float_t sres = 0.25;
247 Float_t res = sd2+sres;
248 return res;
249}
250
251
252Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){
253 //sigma z2 = in digits - angle estimated supposing vertex constraint
753797ce 254 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 255 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
a1ec4d07 256 Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
1c53abe2 257 angular*=angular;
258 angular/=12.;
259 Float_t sres = fParam->GetZSigma()/fZWidth;
260 sres *=sres;
261 Float_t res = angular +sd2+sres;
262 return res;
263}
264
13116aec 265void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
1c53abe2 266AliTPCclusterMI &c)
267{
97f127bb 268 //
269 // k - Make cluster at position k
270 // bins - 2 D array of signals mapped to 1 dimensional array -
271 // max - the number of time bins er one dimension
272 // c - refernce to cluster to be filled
273 //
1c53abe2 274 Int_t i0=k/max; //central pad
275 Int_t j0=k%max; //central time bin
276
277 // set pointers to data
278 //Int_t dummy[5] ={0,0,0,0,0};
13116aec 279 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
1c53abe2 280 for (Int_t di=-2;di<=2;di++){
281 matrix[di+2] = &bins[k+di*max];
1c53abe2 282 }
283 //build matrix with virtual charge
284 Float_t sigmay2= GetSigmaY2(j0);
285 Float_t sigmaz2= GetSigmaZ2(j0);
286
287 Float_t vmatrix[5][5];
288 vmatrix[2][2] = matrix[2][0];
289 c.SetType(0);
6c024a0e 290 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
1c53abe2 291 for (Int_t di =-1;di <=1;di++)
292 for (Int_t dj =-1;dj <=1;dj++){
293 Float_t amp = matrix[di+2][dj];
294 if ( (amp<2) && (fLoop<2)){
295 // if under threshold - calculate virtual charge
296 Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
297 amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
298 if (amp>2) amp = 2;
299 vmatrix[2+di][2+dj]=amp;
300 vmatrix[2+2*di][2+2*dj]=0;
301 if ( (di*dj)!=0){
302 //DIAGONAL ELEMENTS
303 vmatrix[2+2*di][2+dj] =0;
304 vmatrix[2+di][2+2*dj] =0;
305 }
306 continue;
307 }
308 if (amp<4){
309 //if small amplitude - below 2 x threshold - don't consider other one
310 vmatrix[2+di][2+dj]=amp;
311 vmatrix[2+2*di][2+2*dj]=0; // don't take to the account next bin
312 if ( (di*dj)!=0){
313 //DIAGONAL ELEMENTS
314 vmatrix[2+2*di][2+dj] =0;
315 vmatrix[2+di][2+2*dj] =0;
316 }
317 continue;
318 }
319 //if bigger then take everything
320 vmatrix[2+di][2+dj]=amp;
321 vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
322 if ( (di*dj)!=0){
323 //DIAGONAL ELEMENTS
324 vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
325 vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
326 }
327 }
328
329
330
331 Float_t sumw=0;
332 Float_t sumiw=0;
333 Float_t sumi2w=0;
334 Float_t sumjw=0;
335 Float_t sumj2w=0;
336 //
337 for (Int_t i=-2;i<=2;i++)
338 for (Int_t j=-2;j<=2;j++){
339 Float_t amp = vmatrix[i+2][j+2];
340
341 sumw += amp;
342 sumiw += i*amp;
343 sumi2w += i*i*amp;
344 sumjw += j*amp;
345 sumj2w += j*j*amp;
346 }
347 //
348 Float_t meani = sumiw/sumw;
349 Float_t mi2 = sumi2w/sumw-meani*meani;
350 Float_t meanj = sumjw/sumw;
351 Float_t mj2 = sumj2w/sumw-meanj*meanj;
352 //
353 Float_t ry = mi2/sigmay2;
354 Float_t rz = mj2/sigmaz2;
355
356 //
357 if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
194b0609 358 if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) {
359 //
360 //if cluster looks like expected or Unfolding not switched on
361 //standard COG is used
1c53abe2 362 //+1.2 deviation from expected sigma accepted
363 // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
364
365 meani +=i0;
366 meanj +=j0;
367 //set cluster parameters
368 c.SetQ(sumw);
022ee144 369 c.SetPad(meani-2.5);
5bd7ed29 370 c.SetTimeBin(meanj-3);
1c53abe2 371 c.SetSigmaY2(mi2);
372 c.SetSigmaZ2(mj2);
022ee144 373 c.SetType(0);
03fe1804 374 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 375 return;
376 }
377 //
378 //unfolding when neccessary
379 //
380
13116aec 381 Float_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3
382 Float_t dummy[7]={0,0,0,0,0,0};
1c53abe2 383 for (Int_t di=-3;di<=3;di++){
384 matrix2[di+3] = &bins[k+di*max];
385 if ((k+di*max)<3) matrix2[di+3] = &dummy[3];
386 if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3];
387 }
388 Float_t vmatrix2[5][5];
389 Float_t sumu;
390 Float_t overlap;
391 UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
392 //
393 // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
394 meani +=i0;
395 meanj +=j0;
396 //set cluster parameters
397 c.SetQ(sumu);
022ee144 398 c.SetPad(meani-2.5);
399 c.SetTimeBin(meanj-3);
1c53abe2 400 c.SetSigmaY2(mi2);
401 c.SetSigmaZ2(mj2);
402 c.SetType(Char_t(overlap)+1);
03fe1804 403 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 404
405 //unfolding 2
406 meani-=i0;
407 meanj-=j0;
1c53abe2 408}
409
410
411
13116aec 412void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
1c53abe2 413 Float_t & sumu, Float_t & overlap )
414{
415 //
416 //unfold cluster from input matrix
417 //data corresponding to cluster writen in recmatrix
418 //output meani and meanj
419
420 //take separatelly y and z
421
422 Float_t sum3i[7] = {0,0,0,0,0,0,0};
423 Float_t sum3j[7] = {0,0,0,0,0,0,0};
424
425 for (Int_t k =0;k<7;k++)
426 for (Int_t l = -1; l<=1;l++){
427 sum3i[k]+=matrix2[k][l];
428 sum3j[k]+=matrix2[l+3][k-3];
429 }
430 Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
431 //
432 //unfold y
433 Float_t sum3wi = 0; //charge minus overlap
434 Float_t sum3wio = 0; //full charge
435 Float_t sum3iw = 0; //sum for mean value
436 for (Int_t dk=-1;dk<=1;dk++){
437 sum3wio+=sum3i[dk+3];
438 if (dk==0){
439 sum3wi+=sum3i[dk+3];
440 }
441 else{
442 Float_t ratio =1;
443 if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
444 sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
445 Float_t xm2 = sum3i[-dk+3];
446 Float_t xm1 = sum3i[+3];
447 Float_t x1 = sum3i[2*dk+3];
448 Float_t x2 = sum3i[3*dk+3];
cc5e9db0 449 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
450 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 451 ratio = w11/(w11+w12);
452 for (Int_t dl=-1;dl<=1;dl++)
453 mratio[dk+1][dl+1] *= ratio;
454 }
455 Float_t amp = sum3i[dk+3]*ratio;
456 sum3wi+=amp;
457 sum3iw+= dk*amp;
458 }
459 }
460 meani = sum3iw/sum3wi;
461 Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
462
463
464
465 //unfold z
466 Float_t sum3wj = 0; //charge minus overlap
467 Float_t sum3wjo = 0; //full charge
468 Float_t sum3jw = 0; //sum for mean value
469 for (Int_t dk=-1;dk<=1;dk++){
470 sum3wjo+=sum3j[dk+3];
471 if (dk==0){
472 sum3wj+=sum3j[dk+3];
473 }
474 else{
475 Float_t ratio =1;
476 if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
477 (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
478 Float_t xm2 = sum3j[-dk+3];
479 Float_t xm1 = sum3j[+3];
480 Float_t x1 = sum3j[2*dk+3];
481 Float_t x2 = sum3j[3*dk+3];
cc5e9db0 482 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
483 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 484 ratio = w11/(w11+w12);
485 for (Int_t dl=-1;dl<=1;dl++)
486 mratio[dl+1][dk+1] *= ratio;
487 }
488 Float_t amp = sum3j[dk+3]*ratio;
489 sum3wj+=amp;
490 sum3jw+= dk*amp;
491 }
492 }
493 meanj = sum3jw/sum3wj;
494 Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
495 overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
496 sumu = (sum3wj+sum3wi)/2.;
497
498 if (overlap ==3) {
499 //if not overlap detected remove everything
500 for (Int_t di =-2; di<=2;di++)
501 for (Int_t dj =-2; dj<=2;dj++){
502 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
503 }
504 }
505 else{
506 for (Int_t di =-1; di<=1;di++)
507 for (Int_t dj =-1; dj<=1;dj++){
508 Float_t ratio =1;
509 if (mratio[di+1][dj+1]==1){
510 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
511 if (TMath::Abs(di)+TMath::Abs(dj)>1){
512 recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
513 recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
514 }
515 recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
516 }
517 else
518 {
519 //if we have overlap in direction
520 recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
521 if (TMath::Abs(di)+TMath::Abs(dj)>1){
cc5e9db0 522 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
1c53abe2 523 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
524 //
cc5e9db0 525 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
1c53abe2 526 recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
527 }
528 else{
529 ratio = recmatrix[di+2][dj+2]/matrix2[3][0];
530 recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
531 }
532 }
533 }
534 }
1c53abe2 535
536}
537
538Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
539{
540 //
541 // estimate max
542 Float_t sumteor= 0;
543 Float_t sumamp = 0;
544
545 for (Int_t di = -1;di<=1;di++)
546 for (Int_t dj = -1;dj<=1;dj++){
547 if (vmatrix[2+di][2+dj]>2){
548 Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
549 sumteor += teor*vmatrix[2+di][2+dj];
550 sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
551 }
552 }
553 Float_t max = sumamp/sumteor;
554 return max;
555}
556
03fe1804 557void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
1c53abe2 558 //
1c53abe2 559 //
022ee144 560 // Transform cluster to the rotated global coordinata
561 // Assign labels to the cluster
562 // add the cluster to the array
563 // for more details - See AliTPCTranform::Transform(x,i,0,1)
564 Float_t meani = c.GetPad();
565 Float_t meanj = c.GetTimeBin();
1c53abe2 566
022ee144 567 Int_t ki = TMath::Nint(meani);
1c53abe2 568 if (ki<0) ki=0;
569 if (ki>=fMaxPad) ki = fMaxPad-1;
022ee144 570 Int_t kj = TMath::Nint(meanj);
1c53abe2 571 if (kj<0) kj=0;
572 if (kj>=fMaxTime-3) kj=fMaxTime-4;
022ee144 573 // ki and kj shifted as integers coordinata
f8aae377 574 if (fRowDig) {
575 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
576 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
577 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
578 }
1c53abe2 579
022ee144 580 c.SetRow(fRow);
581 c.SetDetector(fSector);
1c53abe2 582 Float_t s2 = c.GetSigmaY2();
583 Float_t w=fParam->GetPadPitchWidth(fSector);
1c53abe2 584 c.SetSigmaY2(s2*w*w);
585 s2 = c.GetSigmaZ2();
022ee144 586 c.SetSigmaZ2(s2*fZWidth*fZWidth);
587 //
588 //
589 //
590 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
591 if (!transform) {
592 AliFatal("Tranformations not in calibDB");
593 }
594 Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
595 Int_t i[1]={fSector};
596 transform->Transform(x,i,0,1);
597 c.SetX(x[0]);
598 c.SetY(x[1]);
599 c.SetZ(x[2]);
600 //
601 //
66a2e32f 602 if (fSector%36>17){
603 c.SetY(-c.GetY());
adefcaa6 604 }
508541c7 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
969 (*fDebugStreamer)<<"Digits"<<
970 "sec="<<fSector<<
971 "row="<<iRow<<
972 "pad="<<iPad<<
973 "time="<<iTimeBin<<
974 "sig="<<signal<<
975 "x="<<x[0]<<
976 "y="<<x[1]<<
977 "z="<<x[2]<<
978 "gx="<<gx[0]<<
979 "gy="<<gx[1]<<
980 "gz="<<gx[2]<<
981 "\n";
982 }
983 }
984 }
985 }
986 }
987
22c352f8 988 // Now loop over rows and find clusters
989 for (fRow = 0; fRow < nRows; fRow++) {
22c352f8 990 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
ebe95b38 991 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
22c352f8 992
993 fRx = fParam->GetPadRowRadii(fSector, fRow);
994 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 995 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 996 if (fSector < kNIS)
997 fMaxPad = fParam->GetNPadsLow(fRow);
998 else
999 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 1000 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 1001
1002 fBins = allBins[fRow];
5b33a7f5 1003 fSigBins = allSigBins[fRow];
1004 fNSigBins = allNSigBins[fRow];
22c352f8 1005
940ed1f0 1006 FindClusters(noiseROC);
98ee6d31 1007
ebe95b38 1008 FillRow();
98ee6d31 1009 if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear();
1010 fNclusters += fNcluster;
1011
22c352f8 1012 } // End of loop to find clusters
22c352f8 1013 } // End of loop over sectors
98ee6d31 1014
adefcaa6 1015 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1016 delete [] allBins[iRow];
5b33a7f5 1017 delete [] allSigBins[iRow];
adefcaa6 1018 }
1019 delete [] allBins;
5b33a7f5 1020 delete [] allSigBins;
1021 delete [] allNSigBins;
adefcaa6 1022
9546a7ff 1023 if (rawReader->GetEventId() && fOutput ){
98ee6d31 1024 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
9546a7ff 1025 }
ebe95b38 1026
db3576a7 1027 if(rawReader->GetEventId()) {
98ee6d31 1028 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1029 }
1030
1031 if(fBClonesArray) {
1032 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
db3576a7 1033 }
f8aae377 1034}
1035
940ed1f0 1036void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 1037{
940ed1f0 1038
1039 //
1040 // add virtual charge at the edge
1041 //
7041b196 1042 Double_t kMaxDumpSize = 500000;
ebe95b38 1043 if (!fOutput) {
1044 fBDumpSignal =kFALSE;
1045 }else{
1046 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1047 }
1048
f8aae377 1049 fNcluster=0;
f8aae377 1050 fLoop=1;
a1ec4d07 1051 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 1052 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1053 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1054 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1055 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1056 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1057 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 1058 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1059 Int_t i = fSigBins[iSig];
1060 if (i%fMaxTime<=crtime) continue;
1061 Float_t *b = &fBins[i];
940ed1f0 1062 //absolute custs
03fe1804 1063 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1064 //
1065 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
aba7fdfc 1066 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1067 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 1068 //
1069 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1070 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1071 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1072 //
1073 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1074 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1075 // sigma cuts
1076 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1077 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1078 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1079
d101caf3 1080 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1081 Int_t dummy=0;
1082 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1083
f8aae377 1084 //}
1085 }
8569a2b0 1086}
65c045f0 1087
eea36fac 1088Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1089 //
679026e6 1090 // Currently hack to filter digital noise (15.06.2008)
eea36fac 1091 // To be parameterized in the AliTPCrecoParam
1092 // More inteligent way to be used in future
679026e6 1093 // Acces to the proper pedestal file needed
eea36fac 1094 //
1095 if (cl->GetMax()<400) return kTRUE;
1096 Double_t ratio = cl->GetQ()/cl->GetMax();
679026e6 1097 if (cl->GetMax()>700){
1098 if ((ratio - int(ratio)>0.8)) return kFALSE;
1099 }
eea36fac 1100 if ((ratio - int(ratio)<0.95)) return kTRUE;
1101 return kFALSE;
eea36fac 1102}
1103
65c045f0 1104
940ed1f0 1105Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1106 //
1107 // process signal on given pad - + streaming of additional information in special mode
1108 //
1109 // id[0] - sector
1110 // id[1] - row
adefcaa6 1111 // id[2] - pad
1112
1113 //
1114 // ESTIMATE pedestal and the noise
1115 //
adefcaa6 1116 const Int_t kPedMax = 100;
1117 Float_t max = 0;
1118 Float_t maxPos = 0;
1119 Int_t median = -1;
1120 Int_t count0 = 0;
1121 Int_t count1 = 0;
940ed1f0 1122 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1123 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1124 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 1125 //
adefcaa6 1126 UShort_t histo[kPedMax];
f174362f 1127 //memset(histo,0,kPedMax*sizeof(UShort_t));
b842afde 1128 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
adefcaa6 1129 for (Int_t i=0; i<fMaxTime; i++){
1130 if (signal[i]<=0) continue;
940ed1f0 1131 if (signal[i]>max && i>firstBin) {
65c045f0 1132 max = signal[i];
1133 maxPos = i;
adefcaa6 1134 }
1135 if (signal[i]>kPedMax-1) continue;
1136 histo[int(signal[i]+0.5)]++;
1137 count0++;
65c045f0 1138 }
7fef31a6 1139 //
adefcaa6 1140 for (Int_t i=1; i<kPedMax; i++){
1141 if (count1<count0*0.5) median=i;
1142 count1+=histo[i];
1143 }
1144 // truncated mean
65c045f0 1145 //
7041b196 1146 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1147 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1148 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1149 //
1150 for (Int_t idelta=1; idelta<10; idelta++){
1151 if (median-idelta<=0) continue;
1152 if (median+idelta>kPedMax) continue;
1153 if (count06<0.6*count1){
1154 count06+=histo[median-idelta];
1155 mean06 +=histo[median-idelta]*(median-idelta);
1156 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1157 count06+=histo[median+idelta];
1158 mean06 +=histo[median+idelta]*(median+idelta);
1159 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1160 }
1161 if (count09<0.9*count1){
1162 count09+=histo[median-idelta];
1163 mean09 +=histo[median-idelta]*(median-idelta);
1164 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1165 count09+=histo[median+idelta];
1166 mean09 +=histo[median+idelta]*(median+idelta);
1167 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1168 }
adefcaa6 1169 if (count10<0.95*count1){
1170 count10+=histo[median-idelta];
1171 mean +=histo[median-idelta]*(median-idelta);
1172 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1173 count10+=histo[median+idelta];
1174 mean +=histo[median+idelta]*(median+idelta);
1175 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1176 }
1177 }
3f0af4ba 1178 if (count10) {
1179 mean /=count10;
1180 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1181 }
1182 if (count06) {
1183 mean06/=count06;
1184 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1185 }
1186 if (count09) {
1187 mean09/=count09;
1188 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1189 }
940ed1f0 1190 rmsEvent = rms09;
adefcaa6 1191 //
940ed1f0 1192 pedestalEvent = median;
adefcaa6 1193 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1194 //
1195 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1196 //
1197 // Dump mean signal info
1198 //
1199 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1200 "TimeStamp="<<fTimeStamp<<
1201 "EventType="<<fEventType<<
adefcaa6 1202 "Sector="<<uid[0]<<
1203 "Row="<<uid[1]<<
1204 "Pad="<<uid[2]<<
1205 "Max="<<max<<
1206 "MaxPos="<<maxPos<<
1207 //
1208 "Median="<<median<<
1209 "Mean="<<mean<<
1210 "RMS="<<rms<<
1211 "Mean06="<<mean06<<
1212 "RMS06="<<rms06<<
1213 "Mean09="<<mean09<<
1214 "RMS09="<<rms09<<
940ed1f0 1215 "RMSCalib="<<rmsCalib<<
1216 "PedCalib="<<pedestalCalib<<
adefcaa6 1217 "\n";
1218 //
1219 // fill pedestal histogram
1220 //
65c045f0 1221 //
adefcaa6 1222 //
1223 //
1224 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1225 Float_t *dsignal = new Float_t[nchannels];
1226 Float_t *dtime = new Float_t[nchannels];
1227 for (Int_t i=0; i<nchannels; i++){
1228 dtime[i] = i;
1229 dsignal[i] = signal[i];
1230 }
03fe1804 1231
11441748 1232 TGraph * graph=0;
7fef31a6 1233 //
a525cc34 1234 // Big signals dumping
1235 //
1236 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1237 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1238 "TimeStamp="<<fTimeStamp<<
1239 "EventType="<<fEventType<<
1240 "Sector="<<uid[0]<<
1241 "Row="<<uid[1]<<
1242 "Pad="<<uid[2]<<
1243 "Graph="<<graph<<
1244 "Max="<<max<<
1245 "MaxPos="<<maxPos<<
1246 //
1247 "Median="<<median<<
1248 "Mean="<<mean<<
1249 "RMS="<<rms<<
1250 "Mean06="<<mean06<<
1251 "RMS06="<<rms06<<
1252 "Mean09="<<mean09<<
1253 "RMS09="<<rms09<<
1254 "\n";
1255 delete graph;
7fef31a6 1256
65c045f0 1257 delete [] dsignal;
1258 delete [] dtime;
940ed1f0 1259 if (rms06>fRecoParam->GetMaxNoise()) {
1260 pedestalEvent+=1024.;
1261 return 1024+median; // sign noisy channel in debug mode
1262 }
65c045f0 1263 return median;
1264}
1265
1266
1267
03fe1804 1268
03fe1804 1269