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