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