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