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