]> 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
using non-uniform field option for tracking (Jacek)
[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
756 if(fBClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
22c352f8 757
98ee6d31 758 // reset counter
759 fNclusters = 0;
88cb7938 760
daac2a58 761 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 762 const Int_t kNIS = fParam->GetNInnerSector();
763 const Int_t kNOS = fParam->GetNOuterSector();
764 const Int_t kNS = kNIS + kNOS;
765 fZWidth = fParam->GetZWidth();
766 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 767 //
768 //alocate memory for sector - maximal case
769 //
22c352f8 770 Float_t** allBins = NULL;
5b33a7f5 771 Int_t** allSigBins = NULL;
772 Int_t* allNSigBins = NULL;
adefcaa6 773 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
774 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
775 allBins = new Float_t*[nRowsMax];
5b33a7f5 776 allSigBins = new Int_t*[nRowsMax];
777 allNSigBins = new Int_t[nRowsMax];
adefcaa6 778 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
779 //
780 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
781 allBins[iRow] = new Float_t[maxBin];
adefcaa6 782 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 783 allSigBins[iRow] = new Int_t[maxBin];
784 allNSigBins[iRow]=0;
adefcaa6 785 }
786 //
22c352f8 787 // Loop over sectors
adefcaa6 788 //
22c352f8 789 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 790
22c352f8 791 Int_t nRows = 0;
792 Int_t nDDLs = 0, indexDDL = 0;
793 if (fSector < kNIS) {
794 nRows = fParam->GetNRowLow();
795 fSign = (fSector < kNIS/2) ? 1 : -1;
796 nDDLs = 2;
797 indexDDL = fSector * 2;
798 }
799 else {
800 nRows = fParam->GetNRowUp();
801 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
802 nDDLs = 4;
803 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
804 }
805
98ee6d31 806 // load the raw data for corresponding DDLs
807 rawReader->Reset();
808 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
809
810 // select only good sector
811 input.Next();
812 if(input.GetSector() != fSector) continue;
813
814 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
815 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
816 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
817 //check the presence of the calibration
818 if (!noiseROC ||!pedestalROC ) {
819 AliError(Form("Missing calibration per sector\t%d\n",fSector));
820 continue;
821 }
822
22c352f8 823 for (Int_t iRow = 0; iRow < nRows; iRow++) {
824 Int_t maxPad;
825 if (fSector < kNIS)
826 maxPad = fParam->GetNPadsLow(iRow);
827 else
828 maxPad = fParam->GetNPadsUp(iRow);
829
830 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
22c352f8 831 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 832 allNSigBins[iRow] = 0;
22c352f8 833 }
834
adefcaa6 835 Int_t digCounter=0;
22c352f8 836 // Begin loop over altro data
adefcaa6 837 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
838 Float_t gain =1;
839 Int_t lastPad=-1;
98ee6d31 840
841 input.Reset();
22c352f8 842 while (input.Next()) {
22c352f8 843 if (input.GetSector() != fSector)
844 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
845
22c352f8 846
847 Int_t iRow = input.GetRow();
9546a7ff 848 if (iRow < 0){
849 continue;
850 }
851
3f0af4ba 852 if (iRow < 0 || iRow >= nRows){
853 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 854 iRow, 0, nRows -1));
3f0af4ba 855 continue;
856 }
adefcaa6 857 //pad
858 Int_t iPad = input.GetPad();
3f0af4ba 859 if (iPad < 0 || iPad >= nPadsMax) {
860 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 861 iPad, 0, nPadsMax-1));
3f0af4ba 862 continue;
863 }
adefcaa6 864 if (iPad!=lastPad){
865 gain = gainROC->GetValue(iRow,iPad);
866 lastPad = iPad;
867 }
868 iPad+=3;
869 //time
870 Int_t iTimeBin = input.GetTime();
daac2a58 871 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
872 continue;
22c352f8 873 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 874 iTimeBin, 0, iTimeBin -1));
daac2a58 875 }
adefcaa6 876 iTimeBin+=3;
3f0af4ba 877
adefcaa6 878 //signal
22c352f8 879 Float_t signal = input.GetSignal();
adefcaa6 880 if (!calcPedestal && signal <= zeroSup) continue;
9546a7ff 881
940ed1f0 882 if (!calcPedestal) {
5b33a7f5 883 Int_t bin = iPad*fMaxTime+iTimeBin;
884 allBins[iRow][bin] = signal/gain;
885 allSigBins[iRow][allNSigBins[iRow]++] = bin;
940ed1f0 886 }else{
887 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
888 }
aba7fdfc 889 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 890
891 // Temporary
892 digCounter++;
22c352f8 893 } // End of the loop over altro data
adefcaa6 894 //
895 //
aba7fdfc 896 //
897 //
22c352f8 898 // Now loop over rows and perform pedestal subtraction
adefcaa6 899 if (digCounter==0) continue;
aba7fdfc 900 // if (calcPedestal) {
901 if (kTRUE) {
22c352f8 902 for (Int_t iRow = 0; iRow < nRows; iRow++) {
903 Int_t maxPad;
904 if (fSector < kNIS)
905 maxPad = fParam->GetNPadsLow(iRow);
906 else
907 maxPad = fParam->GetNPadsUp(iRow);
908
940ed1f0 909 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
aba7fdfc 910 //
911 // Temporary fix for data production - !!!! MARIAN
912 // The noise calibration should take mean and RMS - currently the Gaussian fit used
913 // In case of double peak - the pad should be rejected
914 //
915 // Line mean - if more than given digits over threshold - make a noise calculation
916 // and pedestal substration
917 if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
918 //
940ed1f0 919 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 920 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 921 //Float_t pedestal = TMath::Median(fMaxTime, p);
922 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 923 // calib values
924 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
925 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
926 Double_t rmsEvent = rmsCalib;
927 Double_t pedestalEvent = pedestalCalib;
928 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
929 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
930 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
931
932 //
22c352f8 933 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
5b33a7f5 934 Int_t bin = iPad*fMaxTime+iTimeBin;
935 allBins[iRow][bin] -= pedestalEvent;
7fef31a6 936 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
5b33a7f5 937 allBins[iRow][bin] = 0;
7fef31a6 938 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
5b33a7f5 939 allBins[iRow][bin] = 0;
22c352f8 940 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
5b33a7f5 941 allBins[iRow][bin] = 0;
942 if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
943 allBins[iRow][bin] = 0;
944 if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
22c352f8 945 }
946 }
f8aae377 947 }
22c352f8 948 }
eea36fac 949
950 if (AliTPCReconstructor::StreamLevel()>3) {
951 for (Int_t iRow = 0; iRow < nRows; iRow++) {
952 Int_t maxPad;
953 if (fSector < kNIS)
954 maxPad = fParam->GetNPadsLow(iRow);
955 else
956 maxPad = fParam->GetNPadsUp(iRow);
957
958 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
959 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
960 Int_t bin = iPad*fMaxTime+iTimeBin;
961 Float_t signal = allBins[iRow][bin];
962 if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
963 Double_t x[]={iRow,iPad-3,iTimeBin-3};
964 Int_t i[]={fSector};
965 AliTPCTransform trafo;
966 trafo.Transform(x,i,0,1);
967 Double_t gx[3]={x[0],x[1],x[2]};
968 trafo.RotatedGlobal2Global(fSector,gx);
969
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";
983 }
984 }
985 }
986 }
987 }
988
22c352f8 989 // Now loop over rows and find clusters
990 for (fRow = 0; fRow < nRows; fRow++) {
22c352f8 991 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
ebe95b38 992 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
22c352f8 993
994 fRx = fParam->GetPadRowRadii(fSector, fRow);
995 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 996 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 997 if (fSector < kNIS)
998 fMaxPad = fParam->GetNPadsLow(fRow);
999 else
1000 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 1001 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 1002
1003 fBins = allBins[fRow];
5b33a7f5 1004 fSigBins = allSigBins[fRow];
1005 fNSigBins = allNSigBins[fRow];
22c352f8 1006
940ed1f0 1007 FindClusters(noiseROC);
98ee6d31 1008
ebe95b38 1009 FillRow();
98ee6d31 1010 if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear();
1011 fNclusters += fNcluster;
1012
22c352f8 1013 } // End of loop to find clusters
22c352f8 1014 } // End of loop over sectors
98ee6d31 1015
adefcaa6 1016 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1017 delete [] allBins[iRow];
5b33a7f5 1018 delete [] allSigBins[iRow];
adefcaa6 1019 }
1020 delete [] allBins;
5b33a7f5 1021 delete [] allSigBins;
1022 delete [] allNSigBins;
adefcaa6 1023
9546a7ff 1024 if (rawReader->GetEventId() && fOutput ){
98ee6d31 1025 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
9546a7ff 1026 }
ebe95b38 1027
db3576a7 1028 if(rawReader->GetEventId()) {
98ee6d31 1029 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1030 }
1031
1032 if(fBClonesArray) {
1033 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
db3576a7 1034 }
f8aae377 1035}
1036
940ed1f0 1037void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 1038{
940ed1f0 1039
1040 //
1041 // add virtual charge at the edge
1042 //
7041b196 1043 Double_t kMaxDumpSize = 500000;
ebe95b38 1044 if (!fOutput) {
1045 fBDumpSignal =kFALSE;
1046 }else{
1047 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1048 }
1049
f8aae377 1050 fNcluster=0;
f8aae377 1051 fLoop=1;
a1ec4d07 1052 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 1053 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1054 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1055 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1056 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1057 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1058 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 1059 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1060 Int_t i = fSigBins[iSig];
1061 if (i%fMaxTime<=crtime) continue;
1062 Float_t *b = &fBins[i];
940ed1f0 1063 //absolute custs
03fe1804 1064 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1065 //
1066 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
aba7fdfc 1067 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1068 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 1069 //
1070 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1071 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1072 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1073 //
1074 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1075 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1076 // sigma cuts
1077 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1078 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1079 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1080
d101caf3 1081 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1082 Int_t dummy=0;
1083 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1084
f8aae377 1085 //}
1086 }
8569a2b0 1087}
65c045f0 1088
eea36fac 1089Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1090 //
679026e6 1091 // Currently hack to filter digital noise (15.06.2008)
eea36fac 1092 // To be parameterized in the AliTPCrecoParam
1093 // More inteligent way to be used in future
679026e6 1094 // Acces to the proper pedestal file needed
eea36fac 1095 //
1096 if (cl->GetMax()<400) return kTRUE;
1097 Double_t ratio = cl->GetQ()/cl->GetMax();
679026e6 1098 if (cl->GetMax()>700){
1099 if ((ratio - int(ratio)>0.8)) return kFALSE;
1100 }
eea36fac 1101 if ((ratio - int(ratio)<0.95)) return kTRUE;
1102 return kFALSE;
eea36fac 1103}
1104
65c045f0 1105
940ed1f0 1106Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1107 //
1108 // process signal on given pad - + streaming of additional information in special mode
1109 //
1110 // id[0] - sector
1111 // id[1] - row
adefcaa6 1112 // id[2] - pad
1113
1114 //
1115 // ESTIMATE pedestal and the noise
1116 //
adefcaa6 1117 const Int_t kPedMax = 100;
1118 Float_t max = 0;
1119 Float_t maxPos = 0;
1120 Int_t median = -1;
1121 Int_t count0 = 0;
1122 Int_t count1 = 0;
940ed1f0 1123 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1124 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1125 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 1126 //
adefcaa6 1127 UShort_t histo[kPedMax];
f174362f 1128 //memset(histo,0,kPedMax*sizeof(UShort_t));
b842afde 1129 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
adefcaa6 1130 for (Int_t i=0; i<fMaxTime; i++){
1131 if (signal[i]<=0) continue;
940ed1f0 1132 if (signal[i]>max && i>firstBin) {
65c045f0 1133 max = signal[i];
1134 maxPos = i;
adefcaa6 1135 }
1136 if (signal[i]>kPedMax-1) continue;
1137 histo[int(signal[i]+0.5)]++;
1138 count0++;
65c045f0 1139 }
7fef31a6 1140 //
adefcaa6 1141 for (Int_t i=1; i<kPedMax; i++){
1142 if (count1<count0*0.5) median=i;
1143 count1+=histo[i];
1144 }
1145 // truncated mean
65c045f0 1146 //
7041b196 1147 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1148 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1149 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1150 //
1151 for (Int_t idelta=1; idelta<10; idelta++){
1152 if (median-idelta<=0) continue;
1153 if (median+idelta>kPedMax) continue;
1154 if (count06<0.6*count1){
1155 count06+=histo[median-idelta];
1156 mean06 +=histo[median-idelta]*(median-idelta);
1157 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1158 count06+=histo[median+idelta];
1159 mean06 +=histo[median+idelta]*(median+idelta);
1160 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1161 }
1162 if (count09<0.9*count1){
1163 count09+=histo[median-idelta];
1164 mean09 +=histo[median-idelta]*(median-idelta);
1165 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1166 count09+=histo[median+idelta];
1167 mean09 +=histo[median+idelta]*(median+idelta);
1168 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1169 }
adefcaa6 1170 if (count10<0.95*count1){
1171 count10+=histo[median-idelta];
1172 mean +=histo[median-idelta]*(median-idelta);
1173 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1174 count10+=histo[median+idelta];
1175 mean +=histo[median+idelta]*(median+idelta);
1176 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1177 }
1178 }
3f0af4ba 1179 if (count10) {
1180 mean /=count10;
1181 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1182 }
1183 if (count06) {
1184 mean06/=count06;
1185 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1186 }
1187 if (count09) {
1188 mean09/=count09;
1189 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1190 }
940ed1f0 1191 rmsEvent = rms09;
adefcaa6 1192 //
940ed1f0 1193 pedestalEvent = median;
adefcaa6 1194 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1195 //
1196 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1197 //
1198 // Dump mean signal info
1199 //
1200 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1201 "TimeStamp="<<fTimeStamp<<
1202 "EventType="<<fEventType<<
adefcaa6 1203 "Sector="<<uid[0]<<
1204 "Row="<<uid[1]<<
1205 "Pad="<<uid[2]<<
1206 "Max="<<max<<
1207 "MaxPos="<<maxPos<<
1208 //
1209 "Median="<<median<<
1210 "Mean="<<mean<<
1211 "RMS="<<rms<<
1212 "Mean06="<<mean06<<
1213 "RMS06="<<rms06<<
1214 "Mean09="<<mean09<<
1215 "RMS09="<<rms09<<
940ed1f0 1216 "RMSCalib="<<rmsCalib<<
1217 "PedCalib="<<pedestalCalib<<
adefcaa6 1218 "\n";
1219 //
1220 // fill pedestal histogram
1221 //
65c045f0 1222 //
adefcaa6 1223 //
1224 //
1225 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1226 Float_t *dsignal = new Float_t[nchannels];
1227 Float_t *dtime = new Float_t[nchannels];
1228 for (Int_t i=0; i<nchannels; i++){
1229 dtime[i] = i;
1230 dsignal[i] = signal[i];
1231 }
03fe1804 1232
11441748 1233 TGraph * graph=0;
7fef31a6 1234 //
a525cc34 1235 // Big signals dumping
1236 //
1237 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1238 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1239 "TimeStamp="<<fTimeStamp<<
1240 "EventType="<<fEventType<<
1241 "Sector="<<uid[0]<<
1242 "Row="<<uid[1]<<
1243 "Pad="<<uid[2]<<
1244 "Graph="<<graph<<
1245 "Max="<<max<<
1246 "MaxPos="<<maxPos<<
1247 //
1248 "Median="<<median<<
1249 "Mean="<<mean<<
1250 "RMS="<<rms<<
1251 "Mean06="<<mean06<<
1252 "RMS06="<<rms06<<
1253 "Mean09="<<mean09<<
1254 "RMS09="<<rms09<<
1255 "\n";
1256 delete graph;
7fef31a6 1257
65c045f0 1258 delete [] dsignal;
1259 delete [] dtime;
940ed1f0 1260 if (rms06>fRecoParam->GetMaxNoise()) {
1261 pedestalEvent+=1024.;
1262 return 1024+median; // sign noisy channel in debug mode
1263 }
65c045f0 1264 return median;
1265}
1266
1267
1268
03fe1804 1269
03fe1804 1270