]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCclustererMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bugfix (Matthias)
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
CommitLineData
1c53abe2 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
88cb7938 16/* $Id$ */
1c53abe2 17
18//-------------------------------------------------------
19// Implementation of the TPC clusterer
20//
ebe95b38 21// 1. The Input data for reconstruction - Options
22// 1.a Simulated data - TTree - invoked Digits2Clusters()
23// 1.b Raw data - Digits2Clusters(AliRawReader* rawReader);
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 "AliTPCClustersArray.h"
73#include "AliTPCClustersRow.h"
74#include "AliTPCParam.h"
75#include "AliTPCRawStream.h"
7c9528d8 76#include "AliTPCRawStreamV3.h"
a1e17193 77#include "AliTPCRecoParam.h"
78#include "AliTPCReconstructor.h"
79#include "AliTPCcalibDB.h"
80#include "AliTPCclusterInfo.h"
81#include "AliTPCclusterMI.h"
022ee144 82#include "AliTPCTransform.h"
a1e17193 83#include "AliTPCclustererMI.h"
13116aec 84
1c53abe2 85ClassImp(AliTPCclustererMI)
86
87
88
2f32844a 89AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
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//______________________________________________________________
65c045f0 171AliTPCclustererMI::~AliTPCclustererMI(){
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
1c53abe2 202void AliTPCclustererMI::SetInput(TTree * tree)
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
215void AliTPCclustererMI::SetOutput(TTree * tree)
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
ebe95b38 229void AliTPCclustererMI::FillRow(){
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
1c53abe2 244Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
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
255Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){
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
13116aec 268void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
1c53abe2 269AliTPCclusterMI &c)
270{
97f127bb 271 //
272 // k - Make cluster at position k
273 // bins - 2 D array of signals mapped to 1 dimensional array -
274 // max - the number of time bins er one dimension
275 // c - refernce to cluster to be filled
276 //
1c53abe2 277 Int_t i0=k/max; //central pad
278 Int_t j0=k%max; //central time bin
279
280 // set pointers to data
281 //Int_t dummy[5] ={0,0,0,0,0};
13116aec 282 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
1c53abe2 283 for (Int_t di=-2;di<=2;di++){
284 matrix[di+2] = &bins[k+di*max];
1c53abe2 285 }
286 //build matrix with virtual charge
287 Float_t sigmay2= GetSigmaY2(j0);
288 Float_t sigmaz2= GetSigmaZ2(j0);
289
290 Float_t vmatrix[5][5];
291 vmatrix[2][2] = matrix[2][0];
292 c.SetType(0);
6c024a0e 293 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
1c53abe2 294 for (Int_t di =-1;di <=1;di++)
295 for (Int_t dj =-1;dj <=1;dj++){
296 Float_t amp = matrix[di+2][dj];
297 if ( (amp<2) && (fLoop<2)){
298 // if under threshold - calculate virtual charge
299 Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
300 amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
301 if (amp>2) amp = 2;
302 vmatrix[2+di][2+dj]=amp;
303 vmatrix[2+2*di][2+2*dj]=0;
304 if ( (di*dj)!=0){
305 //DIAGONAL ELEMENTS
306 vmatrix[2+2*di][2+dj] =0;
307 vmatrix[2+di][2+2*dj] =0;
308 }
309 continue;
310 }
311 if (amp<4){
312 //if small amplitude - below 2 x threshold - don't consider other one
313 vmatrix[2+di][2+dj]=amp;
314 vmatrix[2+2*di][2+2*dj]=0; // don't take to the account next bin
315 if ( (di*dj)!=0){
316 //DIAGONAL ELEMENTS
317 vmatrix[2+2*di][2+dj] =0;
318 vmatrix[2+di][2+2*dj] =0;
319 }
320 continue;
321 }
322 //if bigger then take everything
323 vmatrix[2+di][2+dj]=amp;
324 vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
325 if ( (di*dj)!=0){
326 //DIAGONAL ELEMENTS
327 vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
328 vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
329 }
330 }
331
332
333
334 Float_t sumw=0;
335 Float_t sumiw=0;
336 Float_t sumi2w=0;
337 Float_t sumjw=0;
338 Float_t sumj2w=0;
339 //
340 for (Int_t i=-2;i<=2;i++)
341 for (Int_t j=-2;j<=2;j++){
342 Float_t amp = vmatrix[i+2][j+2];
343
344 sumw += amp;
345 sumiw += i*amp;
346 sumi2w += i*i*amp;
347 sumjw += j*amp;
348 sumj2w += j*j*amp;
349 }
350 //
351 Float_t meani = sumiw/sumw;
352 Float_t mi2 = sumi2w/sumw-meani*meani;
353 Float_t meanj = sumjw/sumw;
354 Float_t mj2 = sumj2w/sumw-meanj*meanj;
355 //
356 Float_t ry = mi2/sigmay2;
357 Float_t rz = mj2/sigmaz2;
358
359 //
360 if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
94bf739c 361 if ( ((ry <1.2) && (rz<1.2)) || (!fRecoParam->GetDoUnfold())) {
194b0609 362 //
363 //if cluster looks like expected or Unfolding not switched on
364 //standard COG is used
1c53abe2 365 //+1.2 deviation from expected sigma accepted
366 // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
367
368 meani +=i0;
369 meanj +=j0;
370 //set cluster parameters
371 c.SetQ(sumw);
022ee144 372 c.SetPad(meani-2.5);
5bd7ed29 373 c.SetTimeBin(meanj-3);
1c53abe2 374 c.SetSigmaY2(mi2);
375 c.SetSigmaZ2(mj2);
022ee144 376 c.SetType(0);
03fe1804 377 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 378 return;
379 }
380 //
381 //unfolding when neccessary
382 //
383
13116aec 384 Float_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3
385 Float_t dummy[7]={0,0,0,0,0,0};
1c53abe2 386 for (Int_t di=-3;di<=3;di++){
387 matrix2[di+3] = &bins[k+di*max];
388 if ((k+di*max)<3) matrix2[di+3] = &dummy[3];
389 if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3];
390 }
391 Float_t vmatrix2[5][5];
392 Float_t sumu;
393 Float_t overlap;
394 UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
395 //
396 // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
397 meani +=i0;
398 meanj +=j0;
399 //set cluster parameters
400 c.SetQ(sumu);
022ee144 401 c.SetPad(meani-2.5);
402 c.SetTimeBin(meanj-3);
1c53abe2 403 c.SetSigmaY2(mi2);
404 c.SetSigmaZ2(mj2);
405 c.SetType(Char_t(overlap)+1);
03fe1804 406 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 407
408 //unfolding 2
409 meani-=i0;
410 meanj-=j0;
1c53abe2 411}
412
413
414
13116aec 415void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
1c53abe2 416 Float_t & sumu, Float_t & overlap )
417{
418 //
419 //unfold cluster from input matrix
420 //data corresponding to cluster writen in recmatrix
421 //output meani and meanj
422
423 //take separatelly y and z
424
425 Float_t sum3i[7] = {0,0,0,0,0,0,0};
426 Float_t sum3j[7] = {0,0,0,0,0,0,0};
427
428 for (Int_t k =0;k<7;k++)
429 for (Int_t l = -1; l<=1;l++){
430 sum3i[k]+=matrix2[k][l];
431 sum3j[k]+=matrix2[l+3][k-3];
432 }
433 Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
434 //
435 //unfold y
436 Float_t sum3wi = 0; //charge minus overlap
437 Float_t sum3wio = 0; //full charge
438 Float_t sum3iw = 0; //sum for mean value
439 for (Int_t dk=-1;dk<=1;dk++){
440 sum3wio+=sum3i[dk+3];
441 if (dk==0){
442 sum3wi+=sum3i[dk+3];
443 }
444 else{
445 Float_t ratio =1;
446 if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
94bf739c 447 (sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 )){
1c53abe2 448 Float_t xm2 = sum3i[-dk+3];
449 Float_t xm1 = sum3i[+3];
450 Float_t x1 = sum3i[2*dk+3];
451 Float_t x2 = sum3i[3*dk+3];
cc5e9db0 452 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
453 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 454 ratio = w11/(w11+w12);
455 for (Int_t dl=-1;dl<=1;dl++)
456 mratio[dk+1][dl+1] *= ratio;
457 }
458 Float_t amp = sum3i[dk+3]*ratio;
459 sum3wi+=amp;
460 sum3iw+= dk*amp;
461 }
462 }
463 meani = sum3iw/sum3wi;
464 Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
465
466
467
468 //unfold z
469 Float_t sum3wj = 0; //charge minus overlap
470 Float_t sum3wjo = 0; //full charge
471 Float_t sum3jw = 0; //sum for mean value
472 for (Int_t dk=-1;dk<=1;dk++){
473 sum3wjo+=sum3j[dk+3];
474 if (dk==0){
475 sum3wj+=sum3j[dk+3];
476 }
477 else{
478 Float_t ratio =1;
479 if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
480 (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
481 Float_t xm2 = sum3j[-dk+3];
482 Float_t xm1 = sum3j[+3];
483 Float_t x1 = sum3j[2*dk+3];
484 Float_t x2 = sum3j[3*dk+3];
cc5e9db0 485 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
486 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 487 ratio = w11/(w11+w12);
488 for (Int_t dl=-1;dl<=1;dl++)
489 mratio[dl+1][dk+1] *= ratio;
490 }
491 Float_t amp = sum3j[dk+3]*ratio;
492 sum3wj+=amp;
493 sum3jw+= dk*amp;
494 }
495 }
496 meanj = sum3jw/sum3wj;
497 Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
498 overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
499 sumu = (sum3wj+sum3wi)/2.;
500
501 if (overlap ==3) {
502 //if not overlap detected remove everything
503 for (Int_t di =-2; di<=2;di++)
504 for (Int_t dj =-2; dj<=2;dj++){
505 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
506 }
507 }
508 else{
509 for (Int_t di =-1; di<=1;di++)
510 for (Int_t dj =-1; dj<=1;dj++){
511 Float_t ratio =1;
512 if (mratio[di+1][dj+1]==1){
513 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
514 if (TMath::Abs(di)+TMath::Abs(dj)>1){
515 recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
516 recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
517 }
518 recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
519 }
520 else
521 {
522 //if we have overlap in direction
523 recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
524 if (TMath::Abs(di)+TMath::Abs(dj)>1){
cc5e9db0 525 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
1c53abe2 526 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
527 //
cc5e9db0 528 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
1c53abe2 529 recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
530 }
531 else{
532 ratio = recmatrix[di+2][dj+2]/matrix2[3][0];
533 recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
534 }
535 }
536 }
537 }
1c53abe2 538
539}
540
541Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
542{
543 //
544 // estimate max
545 Float_t sumteor= 0;
546 Float_t sumamp = 0;
547
548 for (Int_t di = -1;di<=1;di++)
549 for (Int_t dj = -1;dj<=1;dj++){
550 if (vmatrix[2+di][2+dj]>2){
551 Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
552 sumteor += teor*vmatrix[2+di][2+dj];
553 sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
554 }
555 }
556 Float_t max = sumamp/sumteor;
557 return max;
558}
559
12ed89e4 560void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * /*matrix*/, Int_t /*pos*/){
1c53abe2 561 //
1c53abe2 562 //
022ee144 563 // Transform cluster to the rotated global coordinata
564 // Assign labels to the cluster
565 // add the cluster to the array
566 // for more details - See AliTPCTranform::Transform(x,i,0,1)
567 Float_t meani = c.GetPad();
568 Float_t meanj = c.GetTimeBin();
1c53abe2 569
022ee144 570 Int_t ki = TMath::Nint(meani);
1c53abe2 571 if (ki<0) ki=0;
572 if (ki>=fMaxPad) ki = fMaxPad-1;
022ee144 573 Int_t kj = TMath::Nint(meanj);
1c53abe2 574 if (kj<0) kj=0;
575 if (kj>=fMaxTime-3) kj=fMaxTime-4;
022ee144 576 // ki and kj shifted as integers coordinata
f8aae377 577 if (fRowDig) {
578 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
579 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
580 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
581 }
1c53abe2 582
022ee144 583 c.SetRow(fRow);
584 c.SetDetector(fSector);
1c53abe2 585 Float_t s2 = c.GetSigmaY2();
586 Float_t w=fParam->GetPadPitchWidth(fSector);
1c53abe2 587 c.SetSigmaY2(s2*w*w);
588 s2 = c.GetSigmaZ2();
022ee144 589 c.SetSigmaZ2(s2*fZWidth*fZWidth);
590 //
591 //
592 //
239f39b1 593
022ee144 594 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
595 if (!transform) {
239f39b1 596 AliFatal("Tranformations not in calibDB");
35e039dc 597 return;
022ee144 598 }
239f39b1 599 transform->SetCurrentRecoParam((AliTPCRecoParam*)fRecoParam);
022ee144 600 Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
601 Int_t i[1]={fSector};
602 transform->Transform(x,i,0,1);
603 c.SetX(x[0]);
604 c.SetY(x[1]);
605 c.SetZ(x[2]);
606 //
607 //
1c53abe2 608 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
1c53abe2 609 c.SetType(-(c.GetType()+3)); //edge clusters
610 }
611 if (fLoop==2) c.SetType(100);
612
98ee6d31 613 // select output
614 TClonesArray * arr = 0;
615 AliTPCclusterMI * cl = 0;
616
617 if(fBClonesArray==kFALSE) {
618 arr = fRowCl->GetArray();
619 cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
620 } else {
621 cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c);
622 }
623
d101caf3 624 // if (fRecoParam->DumpSignal() &&matrix ) {
625// Int_t nbins=0;
626// Float_t *graph =0;
627// if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
628// nbins = fMaxTime;
629// graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
630// }
631// AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
632// cl->SetInfo(info);
633// }
b127a65f 634 if (!fRecoParam->DumpSignal()) {
635 cl->SetInfo(0);
636 }
eea36fac 637
638 if (AliTPCReconstructor::StreamLevel()>1) {
b7a563f9 639 Float_t xyz[3];
640 cl->GetGlobalXYZ(xyz);
eea36fac 641 (*fDebugStreamer)<<"Clusters"<<
642 "Cl.="<<cl<<
b7a563f9 643 "gx="<<xyz[0]<<
644 "gy="<<xyz[1]<<
645 "gz="<<xyz[2]<<
eea36fac 646 "\n";
647 }
1c53abe2 648
649 fNcluster++;
650}
651
652
653//_____________________________________________________________________________
f8aae377 654void AliTPCclustererMI::Digits2Clusters()
1c53abe2 655{
656 //-----------------------------------------------------------------
657 // This is a simple cluster finder.
658 //-----------------------------------------------------------------
1c53abe2 659
f8aae377 660 if (!fInput) {
661 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 662 return;
663 }
b7a563f9 664 fRecoParam = AliTPCReconstructor::GetRecoParam();
665 if (!fRecoParam){
666 AliFatal("Can not get the reconstruction parameters");
667 }
668 if(AliTPCReconstructor::StreamLevel()>5) {
669 AliInfo("Parameter Dumps");
670 fParam->Dump();
671 fRecoParam->Dump();
672 }
a56b3c98 673 fRowDig = NULL;
b7a563f9 674
53a11d77 675 //-----------------------------------------------------------------
676 // Use HLT clusters
677 //-----------------------------------------------------------------
49332743 678 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 679 AliInfo("Using HLT clusters for TPC off-line reconstruction");
680 fZWidth = fParam->GetZWidth();
9abc660b 681 Int_t iResult = ReadHLTClusters();
53a11d77 682
9abc660b 683 // HLT clusters present
bbdbe223 684 if (iResult >= 0 && fNclusters > 0)
685 return;
686
9abc660b 687 // HLT clusters not present
bbdbe223 688 if (iResult < 0 || fNclusters == 0) {
689 if (fUseHLTClusters == 3) {
690 AliError("No HLT clusters present, but requested.");
9abc660b 691 return;
692 }
bbdbe223 693 else {
1598ba75 694 AliInfo("Now trying to read from TPC RAW");
bbdbe223 695 }
9abc660b 696 }
697 // Some other problem during cluster reading
698 else {
699 AliWarning("Some problem while unpacking of HLT clusters.");
700 return;
701 }
49332743 702 } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 703
704 //-----------------------------------------------------------------
705 // Run TPC off-line clusterer
706 //-----------------------------------------------------------------
13116aec 707 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 708 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1c53abe2 709 AliSimDigits digarr, *dummy=&digarr;
710 fRowDig = dummy;
711 fInput->GetBranch("Segment")->SetAddress(&dummy);
712 Stat_t nentries = fInput->GetEntries();
713
daac2a58 714 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 715
1c53abe2 716 Int_t nclusters = 0;
13116aec 717
1c53abe2 718 for (Int_t n=0; n<nentries; n++) {
719 fInput->GetEvent(n);
508541c7 720 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 721 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
722 continue;
723 }
508541c7 724 Int_t row = fRow;
13116aec 725 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 726 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 727 //
1c53abe2 728
ebe95b38 729 fRowCl->SetID(digarr.GetID());
730 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
f8aae377 731 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 732
733
f8aae377 734 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 735 fZWidth = fParam->GetZWidth();
736 if (fSector < kNIS) {
f8aae377 737 fMaxPad = fParam->GetNPadsLow(row);
bbdbe223 738 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 739 fPadLength = fParam->GetPadPitchLength(fSector,row);
740 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 741 } else {
f8aae377 742 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 743 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 744 fPadLength = fParam->GetPadPitchLength(fSector,row);
745 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 746 }
747
748
749 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 750 fBins =new Float_t[fMaxBin];
5b33a7f5 751 fSigBins =new Int_t[fMaxBin];
752 fNSigBins = 0;
13116aec 753 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 754
755 if (digarr.First()) //MI change
756 do {
13116aec 757 Float_t dig=digarr.CurrentDigit();
f8aae377 758 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 759 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 760 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 761 Int_t bin = i*fMaxTime+j;
d6021085 762 if (gain>0){
763 fBins[bin]=dig/gain;
764 }else{
765 fBins[bin]=0;
766 }
5b33a7f5 767 fSigBins[fNSigBins++]=bin;
1c53abe2 768 } while (digarr.Next());
769 digarr.ExpandTrackBuffer();
770
940ed1f0 771 FindClusters(noiseROC);
ebe95b38 772 FillRow();
e7de6656 773 fRowCl->GetArray()->Clear("C");
88cb7938 774 nclusters+=fNcluster;
98ee6d31 775
5b33a7f5 776 delete[] fBins;
777 delete[] fSigBins;
88cb7938 778 }
b7a563f9 779
19dd5b2f 780 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
bbdbe223 781
782 if (fUseHLTClusters == 2 && nclusters == 0) {
783 AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
784
785 fZWidth = fParam->GetZWidth();
786 ReadHLTClusters();
787 }
f8aae377 788}
789
e7034a8b 790void AliTPCclustererMI::ProcessSectorData(){
7c9528d8 791 //
792 // Process the data for the current sector
793 //
794
795 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
796 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
797 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
798 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
799 //check the presence of the calibration
800 if (!noiseROC ||!pedestalROC ) {
801 AliError(Form("Missing calibration per sector\t%d\n",fSector));
802 return;
803 }
804 Int_t nRows=fParam->GetNRow(fSector);
805 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
806 Int_t zeroSup = fParam->GetZeroSup();
807 // if (calcPedestal) {
808 if (kFALSE ) {
809 for (Int_t iRow = 0; iRow < nRows; iRow++) {
810 Int_t maxPad = fParam->GetNPads(fSector, iRow);
811
812 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
813 //
814 // Temporary fix for data production - !!!! MARIAN
815 // The noise calibration should take mean and RMS - currently the Gaussian fit used
816 // In case of double peak - the pad should be rejected
817 //
818 // Line mean - if more than given digits over threshold - make a noise calculation
819 // and pedestal substration
e7034a8b 820 if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue;
7c9528d8 821 //
e7034a8b 822 if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
823 Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3];
7c9528d8 824 //Float_t pedestal = TMath::Median(fMaxTime, p);
825 Int_t id[3] = {fSector, iRow, iPad-3};
826 // calib values
827 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
828 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
829 Double_t rmsEvent = rmsCalib;
830 Double_t pedestalEvent = pedestalCalib;
831 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
832 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
833 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
834
835 //
836 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
837 Int_t bin = iPad*fMaxTime+iTimeBin;
e7034a8b 838 fAllBins[iRow][bin] -= pedestalEvent;
7c9528d8 839 if (iTimeBin < fRecoParam->GetFirstBin())
e7034a8b 840 fAllBins[iRow][bin] = 0;
7c9528d8 841 if (iTimeBin > fRecoParam->GetLastBin())
e7034a8b 842 fAllBins[iRow][bin] = 0;
843 if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
844 fAllBins[iRow][bin] = 0;
845 if (fAllBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
846 fAllBins[iRow][bin] = 0;
847 if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
7c9528d8 848 }
849 }
850 }
851 }
852
6fbe1e5c 853 if (AliTPCReconstructor::StreamLevel()>5) {
7c9528d8 854 for (Int_t iRow = 0; iRow < nRows; iRow++) {
855 Int_t maxPad = fParam->GetNPads(fSector,iRow);
856
857 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
858 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
859 Int_t bin = iPad*fMaxTime+iTimeBin;
e7034a8b 860 Float_t signal = fAllBins[iRow][bin];
7c9528d8 861 if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
862 Double_t x[]={iRow,iPad-3,iTimeBin-3};
863 Int_t i[]={fSector};
864 AliTPCTransform trafo;
865 trafo.Transform(x,i,0,1);
866 Double_t gx[3]={x[0],x[1],x[2]};
867 trafo.RotatedGlobal2Global(fSector,gx);
e7034a8b 868 // fAllSigBins[iRow][fAllNSigBins[iRow]++]
869 Int_t rowsigBins = fAllNSigBins[iRow];
870 Int_t first=fAllSigBins[iRow][0];
7c9528d8 871 Int_t last= 0;
e7034a8b 872 // if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1];
7c9528d8 873
6fbe1e5c 874 if (AliTPCReconstructor::StreamLevel()>5) {
7c9528d8 875 (*fDebugStreamer)<<"Digits"<<
876 "sec="<<fSector<<
877 "row="<<iRow<<
878 "pad="<<iPad<<
879 "time="<<iTimeBin<<
880 "sig="<<signal<<
881 "x="<<x[0]<<
882 "y="<<x[1]<<
883 "z="<<x[2]<<
884 "gx="<<gx[0]<<
885 "gy="<<gx[1]<<
886 "gz="<<gx[2]<<
887 //
888 "rowsigBins="<<rowsigBins<<
889 "first="<<first<<
890 "last="<<last<<
891 "\n";
892 }
893 }
894 }
895 }
896 }
897 }
898
899 // Now loop over rows and find clusters
900 for (fRow = 0; fRow < nRows; fRow++) {
901 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
902 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
903
904 fRx = fParam->GetPadRowRadii(fSector, fRow);
905 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
906 fPadWidth = fParam->GetPadPitchWidth();
907 fMaxPad = fParam->GetNPads(fSector,fRow);
908 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
909
e7034a8b 910 fBins = fAllBins[fRow];
911 fSigBins = fAllSigBins[fRow];
912 fNSigBins = fAllNSigBins[fRow];
7c9528d8 913
914 FindClusters(noiseROC);
915
916 FillRow();
e7de6656 917 if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C");
7c9528d8 918 fNclusters += fNcluster;
919
920 } // End of loop to find clusters
921}
922
923
f8aae377 924void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
925{
926//-----------------------------------------------------------------
22c352f8 927// This is a cluster finder for the TPC raw data.
928// The method assumes NO ordering of the altro channels.
929// The pedestal subtraction can be switched on and off
930// using an option of the TPC reconstructor
f8aae377 931//-----------------------------------------------------------------
b7a563f9 932 fRecoParam = AliTPCReconstructor::GetRecoParam();
933 if (!fRecoParam){
934 AliFatal("Can not get the reconstruction parameters");
935 }
936 if(AliTPCReconstructor::StreamLevel()>5) {
937 AliInfo("Parameter Dumps");
938 fParam->Dump();
939 fRecoParam->Dump();
940 }
f8aae377 941 fRowDig = NULL;
e7034a8b 942
53a11d77 943 //-----------------------------------------------------------------
944 // Use HLT clusters
945 //-----------------------------------------------------------------
49332743 946 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 947 AliInfo("Using HLT clusters for TPC off-line reconstruction");
948 fZWidth = fParam->GetZWidth();
9abc660b 949 Int_t iResult = ReadHLTClusters();
53a11d77 950
9abc660b 951 // HLT clusters present
bbdbe223 952 if (iResult >= 0 && fNclusters > 0)
9abc660b 953 return;
bbdbe223 954
9abc660b 955 // HLT clusters not present
bbdbe223 956 if (iResult < 0 || fNclusters == 0) {
957 if (fUseHLTClusters == 3) {
958 AliError("No HLT clusters present, but requested.");
9abc660b 959 return;
960 }
bbdbe223 961 else {
962 AliInfo("Now trying to read TPC RAW");
963 }
9abc660b 964 }
965 // Some other problem during cluster reading
966 else {
967 AliWarning("Some problem while unpacking of HLT clusters.");
968 return;
969 }
49332743 970 } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 971
972 //-----------------------------------------------------------------
973 // Run TPC off-line clusterer
974 //-----------------------------------------------------------------
22c352f8 975 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
7c9528d8 976 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
977 //
978 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
979 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
980 if (fEventHeader){
981 fTimeStamp = fEventHeader->Get("Timestamp");
982 fEventType = fEventHeader->Get("Type");
fefec90f 983 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
984 transform->SetCurrentTimeStamp(fTimeStamp);
985 transform->SetCurrentRun(rawReader->GetRunNumber());
7c9528d8 986 }
987
988 // creaate one TClonesArray for all clusters
989 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
990 // reset counter
991 fNclusters = 0;
992
993 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
994// const Int_t kNIS = fParam->GetNInnerSector();
995// const Int_t kNOS = fParam->GetNOuterSector();
996// const Int_t kNS = kNIS + kNOS;
997 fZWidth = fParam->GetZWidth();
998 Int_t zeroSup = fParam->GetZeroSup();
999 //
e7034a8b 1000 // Clean-up
7c9528d8 1001 //
e7034a8b 1002 AliTPCROC * roc = AliTPCROC::Instance();
7c9528d8 1003 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1004 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
7c9528d8 1005 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1006 //
1007 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
e7034a8b 1008 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1009 fAllNSigBins[iRow]=0;
7c9528d8 1010 }
1011
1012 Int_t prevSector=-1;
1013 rawReader->Reset();
1014 Int_t digCounter=0;
1015 //
1016 // Loop over DDLs
1017 //
1018 const Int_t kNIS = fParam->GetNInnerSector();
1019 const Int_t kNOS = fParam->GetNOuterSector();
1020 const Int_t kNS = kNIS + kNOS;
1021
1022 for(fSector = 0; fSector < kNS; fSector++) {
1023
1024 Int_t nRows = 0;
1025 Int_t nDDLs = 0, indexDDL = 0;
1026 if (fSector < kNIS) {
1027 nRows = fParam->GetNRowLow();
1028 fSign = (fSector < kNIS/2) ? 1 : -1;
1029 nDDLs = 2;
1030 indexDDL = fSector * 2;
1031 }
1032 else {
1033 nRows = fParam->GetNRowUp();
1034 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1035 nDDLs = 4;
1036 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1037 }
1038
1039 // load the raw data for corresponding DDLs
1040 rawReader->Reset();
1041 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1042
1043 while (input.NextDDL()){
1044 if (input.GetSector() != fSector)
1045 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1046
fefec90f 1047 //Int_t nRows = fParam->GetNRow(fSector);
7c9528d8 1048
1049 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
1050 // Begin loop over altro data
1051 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1052 Float_t gain =1;
1053
1054 //loop over pads
1055 while ( input.NextChannel() ) {
1056 Int_t iRow = input.GetRow();
1057 if (iRow < 0){
1058 continue;
1059 }
1060 if (iRow >= nRows){
1061 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1062 iRow, 0, nRows -1));
1063 continue;
1064 }
1065 //pad
1066 Int_t iPad = input.GetPad();
1067 if (iPad < 0 || iPad >= nPadsMax) {
1068 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1069 iPad, 0, nPadsMax-1));
1070 continue;
1071 }
1072 gain = gainROC->GetValue(iRow,iPad);
1073 iPad+=3;
1074
1075 //loop over bunches
1076 while ( input.NextBunch() ){
1077 Int_t startTbin = (Int_t)input.GetStartTimeBin();
1078 Int_t bunchlength = (Int_t)input.GetBunchLength();
1079 const UShort_t *sig = input.GetSignals();
1080 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1081 Int_t iTimeBin=startTbin-iTime;
1082 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1083 continue;
1084 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1085 iTimeBin, 0, iTimeBin -1));
1086 }
1087 iTimeBin+=3;
1088 //signal
1089 Float_t signal=(Float_t)sig[iTime];
1090 if (!calcPedestal && signal <= zeroSup) continue;
1091
1092 if (!calcPedestal) {
1093 Int_t bin = iPad*fMaxTime+iTimeBin;
1094 if (gain>0){
e7034a8b 1095 fAllBins[iRow][bin] = signal/gain;
7c9528d8 1096 }else{
e7034a8b 1097 fAllBins[iRow][bin] =0;
7c9528d8 1098 }
e7034a8b 1099 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
7c9528d8 1100 }else{
e7034a8b 1101 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
7c9528d8 1102 }
e7034a8b 1103 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
7c9528d8 1104
1105 // Temporary
1106 digCounter++;
1107 }// end loop signals in bunch
1108 }// end loop bunches
1109 } // end loop pads
1110 //
1111 //
1112 //
1113 //
1114 // Now loop over rows and perform pedestal subtraction
1115 if (digCounter==0) continue;
1116 } // End of loop over sectors
1117 //process last sector
1118 if ( digCounter>0 ){
e7034a8b 1119 ProcessSectorData();
7c9528d8 1120 for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1121 Int_t maxPad = fParam->GetNPads(fSector,iRow);
1122 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
e7034a8b 1123 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1124 fAllNSigBins[iRow] = 0;
7c9528d8 1125 }
1126 prevSector=fSector;
1127 digCounter=0;
1128 }
1129 }
1130
7c9528d8 1131 if (rawReader->GetEventId() && fOutput ){
3a4b0aed 1132 Info("Digits2Clusters", "File %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
7c9528d8 1133 }
1134
1135 if(rawReader->GetEventId()) {
3a4b0aed 1136 Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
7c9528d8 1137 }
1138
1139 if(fBClonesArray) {
1140 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1141 }
bbdbe223 1142
1143 if (fUseHLTClusters == 2 && fNclusters == 0) {
1144 AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
1145
1146 fZWidth = fParam->GetZWidth();
1147 ReadHLTClusters();
1148 }
7c9528d8 1149}
1150
1151
1152
1153
1154
1155void AliTPCclustererMI::Digits2ClustersOld
1156(AliRawReader* rawReader)
1157{
1158//-----------------------------------------------------------------
1159// This is a cluster finder for the TPC raw data.
1160// The method assumes NO ordering of the altro channels.
1161// The pedestal subtraction can be switched on and off
1162// using an option of the TPC reconstructor
1163//-----------------------------------------------------------------
1164 fRecoParam = AliTPCReconstructor::GetRecoParam();
1165 if (!fRecoParam){
1166 AliFatal("Can not get the reconstruction parameters");
1167 }
1168 if(AliTPCReconstructor::StreamLevel()>5) {
1169 AliInfo("Parameter Dumps");
1170 fParam->Dump();
1171 fRecoParam->Dump();
1172 }
1173 fRowDig = NULL;
e7034a8b 1174
7c9528d8 1175 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
d6834f5f 1176 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1177 //
1178 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 1179 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1180 if (fEventHeader){
1181 fTimeStamp = fEventHeader->Get("Timestamp");
1182 fEventType = fEventHeader->Get("Type");
1183 }
1184
98ee6d31 1185 // creaate one TClonesArray for all clusters
44adbd4b 1186 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
98ee6d31 1187 // reset counter
1188 fNclusters = 0;
88cb7938 1189
daac2a58 1190 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 1191 const Int_t kNIS = fParam->GetNInnerSector();
1192 const Int_t kNOS = fParam->GetNOuterSector();
1193 const Int_t kNS = kNIS + kNOS;
1194 fZWidth = fParam->GetZWidth();
1195 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 1196 //
e7034a8b 1197 // Clean-up
adefcaa6 1198 //
e7034a8b 1199
1200 AliTPCROC * roc = AliTPCROC::Instance();
adefcaa6 1201 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1202 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
adefcaa6 1203 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1204 //
1205 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
e7034a8b 1206 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1207 fAllNSigBins[iRow]=0;
adefcaa6 1208 }
1209 //
22c352f8 1210 // Loop over sectors
adefcaa6 1211 //
22c352f8 1212 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 1213
22c352f8 1214 Int_t nRows = 0;
1215 Int_t nDDLs = 0, indexDDL = 0;
1216 if (fSector < kNIS) {
1217 nRows = fParam->GetNRowLow();
1218 fSign = (fSector < kNIS/2) ? 1 : -1;
1219 nDDLs = 2;
1220 indexDDL = fSector * 2;
1221 }
1222 else {
1223 nRows = fParam->GetNRowUp();
1224 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1225 nDDLs = 4;
1226 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1227 }
1228
98ee6d31 1229 // load the raw data for corresponding DDLs
1230 rawReader->Reset();
1231 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1232
1233 // select only good sector
35e039dc 1234 if (!input.Next()) continue;
98ee6d31 1235 if(input.GetSector() != fSector) continue;
1236
1237 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
98ee6d31 1238
22c352f8 1239 for (Int_t iRow = 0; iRow < nRows; iRow++) {
1240 Int_t maxPad;
1241 if (fSector < kNIS)
7c9528d8 1242 maxPad = fParam->GetNPadsLow(iRow);
22c352f8 1243 else
7c9528d8 1244 maxPad = fParam->GetNPadsUp(iRow);
22c352f8 1245
1246 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
e7034a8b 1247 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1248 fAllNSigBins[iRow] = 0;
22c352f8 1249 }
1250
adefcaa6 1251 Int_t digCounter=0;
22c352f8 1252 // Begin loop over altro data
adefcaa6 1253 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1254 Float_t gain =1;
1255 Int_t lastPad=-1;
98ee6d31 1256
1257 input.Reset();
22c352f8 1258 while (input.Next()) {
22c352f8 1259 if (input.GetSector() != fSector)
7c9528d8 1260 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1261
22c352f8 1262
22c352f8 1263 Int_t iRow = input.GetRow();
9546a7ff 1264 if (iRow < 0){
7c9528d8 1265 continue;
9546a7ff 1266 }
1267
3f0af4ba 1268 if (iRow < 0 || iRow >= nRows){
7c9528d8 1269 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 1270 iRow, 0, nRows -1));
7c9528d8 1271 continue;
3f0af4ba 1272 }
adefcaa6 1273 //pad
1274 Int_t iPad = input.GetPad();
3f0af4ba 1275 if (iPad < 0 || iPad >= nPadsMax) {
7c9528d8 1276 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 1277 iPad, 0, nPadsMax-1));
7c9528d8 1278 continue;
3f0af4ba 1279 }
adefcaa6 1280 if (iPad!=lastPad){
7c9528d8 1281 gain = gainROC->GetValue(iRow,iPad);
1282 lastPad = iPad;
adefcaa6 1283 }
1284 iPad+=3;
1285 //time
1286 Int_t iTimeBin = input.GetTime();
daac2a58 1287 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
7c9528d8 1288 continue;
1289 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 1290 iTimeBin, 0, iTimeBin -1));
daac2a58 1291 }
adefcaa6 1292 iTimeBin+=3;
3f0af4ba 1293
adefcaa6 1294 //signal
22c352f8 1295 Float_t signal = input.GetSignal();
7c9528d8 1296 if (!calcPedestal && signal <= zeroSup) continue;
9546a7ff 1297
940ed1f0 1298 if (!calcPedestal) {
7c9528d8 1299 Int_t bin = iPad*fMaxTime+iTimeBin;
1300 if (gain>0){
e7034a8b 1301 fAllBins[iRow][bin] = signal/gain;
7c9528d8 1302 }else{
e7034a8b 1303 fAllBins[iRow][bin] =0;
7c9528d8 1304 }
e7034a8b 1305 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
940ed1f0 1306 }else{
e7034a8b 1307 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
940ed1f0 1308 }
e7034a8b 1309 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 1310
1311 // Temporary
1312 digCounter++;
22c352f8 1313 } // End of the loop over altro data
adefcaa6 1314 //
1315 //
aba7fdfc 1316 //
1317 //
22c352f8 1318 // Now loop over rows and perform pedestal subtraction
adefcaa6 1319 if (digCounter==0) continue;
e7034a8b 1320 ProcessSectorData();
22c352f8 1321 } // End of loop over sectors
98ee6d31 1322
9546a7ff 1323 if (rawReader->GetEventId() && fOutput ){
98ee6d31 1324 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
9546a7ff 1325 }
ebe95b38 1326
db3576a7 1327 if(rawReader->GetEventId()) {
98ee6d31 1328 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1329 }
1330
1331 if(fBClonesArray) {
1332 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
db3576a7 1333 }
f8aae377 1334}
1335
940ed1f0 1336void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 1337{
940ed1f0 1338
1339 //
1340 // add virtual charge at the edge
1341 //
7041b196 1342 Double_t kMaxDumpSize = 500000;
ebe95b38 1343 if (!fOutput) {
1344 fBDumpSignal =kFALSE;
1345 }else{
1346 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1347 }
1348
f8aae377 1349 fNcluster=0;
f8aae377 1350 fLoop=1;
a1ec4d07 1351 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 1352 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1353 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1354 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1355 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1356 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1357 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
f47588e0 1358 Int_t useOnePadCluster = fRecoParam->GetUseOnePadCluster();
5b33a7f5 1359 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1360 Int_t i = fSigBins[iSig];
1bfaa9e9 1361 if (i%fMaxTime<=crtime) continue;
5b33a7f5 1362 Float_t *b = &fBins[i];
940ed1f0 1363 //absolute custs
03fe1804 1364 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1365 //
f47588e0 1366 if (useOnePadCluster==0){
1367 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1368 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1369 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1370 }
03fe1804 1371 //
1372 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1373 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1374 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1375 //
1376 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1377 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1378 // sigma cuts
1379 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1380 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1381 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1382
d101caf3 1383 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1384 Int_t dummy=0;
1385 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1386
f8aae377 1387 //}
1388 }
8569a2b0 1389}
65c045f0 1390
eea36fac 1391Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
80efd770 1392 // -- Depricated --
679026e6 1393 // Currently hack to filter digital noise (15.06.2008)
eea36fac 1394 // To be parameterized in the AliTPCrecoParam
1395 // More inteligent way to be used in future
679026e6 1396 // Acces to the proper pedestal file needed
eea36fac 1397 //
1398 if (cl->GetMax()<400) return kTRUE;
1399 Double_t ratio = cl->GetQ()/cl->GetMax();
679026e6 1400 if (cl->GetMax()>700){
1401 if ((ratio - int(ratio)>0.8)) return kFALSE;
1402 }
eea36fac 1403 if ((ratio - int(ratio)<0.95)) return kTRUE;
1404 return kFALSE;
eea36fac 1405}
1406
65c045f0 1407
940ed1f0 1408Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1409 //
1410 // process signal on given pad - + streaming of additional information in special mode
1411 //
1412 // id[0] - sector
1413 // id[1] - row
adefcaa6 1414 // id[2] - pad
1415
1416 //
1417 // ESTIMATE pedestal and the noise
1418 //
adefcaa6 1419 const Int_t kPedMax = 100;
1420 Float_t max = 0;
1421 Float_t maxPos = 0;
1422 Int_t median = -1;
1423 Int_t count0 = 0;
1424 Int_t count1 = 0;
940ed1f0 1425 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1426 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
b7a563f9 1427 Int_t firstBin = fRecoParam->GetFirstBin();
65c045f0 1428 //
adefcaa6 1429 UShort_t histo[kPedMax];
f174362f 1430 //memset(histo,0,kPedMax*sizeof(UShort_t));
b842afde 1431 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
adefcaa6 1432 for (Int_t i=0; i<fMaxTime; i++){
1433 if (signal[i]<=0) continue;
940ed1f0 1434 if (signal[i]>max && i>firstBin) {
65c045f0 1435 max = signal[i];
1436 maxPos = i;
adefcaa6 1437 }
1438 if (signal[i]>kPedMax-1) continue;
1439 histo[int(signal[i]+0.5)]++;
1440 count0++;
65c045f0 1441 }
7fef31a6 1442 //
adefcaa6 1443 for (Int_t i=1; i<kPedMax; i++){
1444 if (count1<count0*0.5) median=i;
1445 count1+=histo[i];
1446 }
1447 // truncated mean
65c045f0 1448 //
7041b196 1449 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1450 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1451 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1452 //
1453 for (Int_t idelta=1; idelta<10; idelta++){
1454 if (median-idelta<=0) continue;
1455 if (median+idelta>kPedMax) continue;
1456 if (count06<0.6*count1){
1457 count06+=histo[median-idelta];
1458 mean06 +=histo[median-idelta]*(median-idelta);
1459 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1460 count06+=histo[median+idelta];
1461 mean06 +=histo[median+idelta]*(median+idelta);
1462 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1463 }
1464 if (count09<0.9*count1){
1465 count09+=histo[median-idelta];
1466 mean09 +=histo[median-idelta]*(median-idelta);
1467 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1468 count09+=histo[median+idelta];
1469 mean09 +=histo[median+idelta]*(median+idelta);
1470 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1471 }
adefcaa6 1472 if (count10<0.95*count1){
1473 count10+=histo[median-idelta];
1474 mean +=histo[median-idelta]*(median-idelta);
1475 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1476 count10+=histo[median+idelta];
1477 mean +=histo[median+idelta]*(median+idelta);
1478 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1479 }
1480 }
3f0af4ba 1481 if (count10) {
1482 mean /=count10;
1483 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1484 }
1485 if (count06) {
1486 mean06/=count06;
1487 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1488 }
1489 if (count09) {
1490 mean09/=count09;
1491 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1492 }
940ed1f0 1493 rmsEvent = rms09;
adefcaa6 1494 //
940ed1f0 1495 pedestalEvent = median;
adefcaa6 1496 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1497 //
1498 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1499 //
1500 // Dump mean signal info
1501 //
b2426624 1502 if (AliTPCReconstructor::StreamLevel()>0) {
1503 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1504 "TimeStamp="<<fTimeStamp<<
1505 "EventType="<<fEventType<<
adefcaa6 1506 "Sector="<<uid[0]<<
1507 "Row="<<uid[1]<<
1508 "Pad="<<uid[2]<<
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<<
940ed1f0 1519 "RMSCalib="<<rmsCalib<<
1520 "PedCalib="<<pedestalCalib<<
adefcaa6 1521 "\n";
b2426624 1522 }
adefcaa6 1523 //
1524 // fill pedestal histogram
1525 //
65c045f0 1526 //
adefcaa6 1527 //
1528 //
1529 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1530 Float_t *dsignal = new Float_t[nchannels];
1531 Float_t *dtime = new Float_t[nchannels];
1532 for (Int_t i=0; i<nchannels; i++){
1533 dtime[i] = i;
1534 dsignal[i] = signal[i];
1535 }
03fe1804 1536
7fef31a6 1537 //
a525cc34 1538 // Big signals dumping
1539 //
b2426624 1540 if (AliTPCReconstructor::StreamLevel()>0) {
b7a563f9 1541 if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin())
a525cc34 1542 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1543 "TimeStamp="<<fTimeStamp<<
1544 "EventType="<<fEventType<<
1545 "Sector="<<uid[0]<<
1546 "Row="<<uid[1]<<
1547 "Pad="<<uid[2]<<
35e039dc 1548 // "Graph="<<graph<<
a525cc34 1549 "Max="<<max<<
1550 "MaxPos="<<maxPos<<
1551 //
1552 "Median="<<median<<
1553 "Mean="<<mean<<
1554 "RMS="<<rms<<
1555 "Mean06="<<mean06<<
1556 "RMS06="<<rms06<<
1557 "Mean09="<<mean09<<
1558 "RMS09="<<rms09<<
1559 "\n";
b2426624 1560 }
7fef31a6 1561
65c045f0 1562 delete [] dsignal;
1563 delete [] dtime;
940ed1f0 1564 if (rms06>fRecoParam->GetMaxNoise()) {
1565 pedestalEvent+=1024.;
1566 return 1024+median; // sign noisy channel in debug mode
1567 }
65c045f0 1568 return median;
1569}
1570
53a11d77 1571Int_t AliTPCclustererMI::ReadHLTClusters()
1572{
1573 //
1574 // read HLT clusters instead of off line custers,
1575 // used in Digits2Clusters
1576 //
1577
6e1c6e9f 1578 if (!fHLTClusterAccess) {
53a11d77 1579 TClass* pCl=NULL;
1580 ROOT::NewFunc_t pNewFunc=NULL;
1581 do {
1582 pCl=TClass::GetClass("AliHLTTPCClusterAccessHLTOUT");
1583 } while (!pCl && gSystem->Load("libAliHLTTPC.so")==0);
1584 if (!pCl || (pNewFunc=pCl->GetNew())==NULL) {
1585 AliError("can not load class description of AliHLTTPCClusterAccessHLTOUT, aborting ...");
1586 return -1;
1587 }
1588
1589 void* p=(*pNewFunc)(NULL);
1590 if (!p) {
1591 AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT");
1592 return -2;
1593 }
6e1c6e9f 1594 fHLTClusterAccess=reinterpret_cast<TObject*>(p);
53a11d77 1595 }
1596
6e1c6e9f 1597 TObject* pClusterAccess=fHLTClusterAccess;
1598
53a11d77 1599 const Int_t kNIS = fParam->GetNInnerSector();
1600 const Int_t kNOS = fParam->GetNOuterSector();
1601 const Int_t kNS = kNIS + kNOS;
1602 fNclusters = 0;
1603
6e1c6e9f 1604 // make sure that all clusters from the previous event are cleared
1605 pClusterAccess->Clear("event");
53a11d77 1606 for(fSector = 0; fSector < kNS; fSector++) {
1607
bbdbe223 1608 Int_t iResult = 1;
53a11d77 1609 TString param("sector="); param+=fSector;
6e1c6e9f 1610 // prepare for next sector
1611 pClusterAccess->Clear("sector");
bbdbe223 1612 pClusterAccess->Execute("read", param, &iResult);
1613 if (iResult < 0) {
1614 return iResult;
1615 AliError("HLT Clusters can not be found");
1616 }
1617
53a11d77 1618 if (pClusterAccess->FindObject("clusterarray")==NULL) {
1619 AliError("HLT clusters requested, but not cluster array not present");
1620 return -4;
1621 }
1622
1623 TClonesArray* clusterArray=dynamic_cast<TClonesArray*>(pClusterAccess->FindObject("clusterarray"));
1624 if (!clusterArray) {
1625 AliError("HLT cluster array is not of class type TClonesArray");
1626 return -5;
1627 }
1628
1629 AliDebug(4,Form("Reading %d clusters from HLT for sector %d", clusterArray->GetEntriesFast(), fSector));
65c045f0 1630
53a11d77 1631 Int_t nClusterSector=0;
1632 Int_t nRows=fParam->GetNRow(fSector);
65c045f0 1633
53a11d77 1634 for (fRow = 0; fRow < nRows; fRow++) {
1635 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
1636 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
1637 fNcluster=0; // reset clusters per row
1638
1639 fRx = fParam->GetPadRowRadii(fSector, fRow);
1640 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
1641 fPadWidth = fParam->GetPadPitchWidth();
1642 fMaxPad = fParam->GetNPads(fSector,fRow);
1643 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
1644
1645 fBins = fAllBins[fRow];
1646 fSigBins = fAllSigBins[fRow];
1647 fNSigBins = fAllNSigBins[fRow];
1648
1649 for (Int_t i=0; i<clusterArray->GetEntriesFast(); i++) {
1650 if (!clusterArray->At(i))
1651 continue;
1652
1653 AliTPCclusterMI* cluster=dynamic_cast<AliTPCclusterMI*>(clusterArray->At(i));
1654 if (!cluster) continue;
1655 if (cluster->GetRow()!=fRow) continue;
1656 nClusterSector++;
1657 AddCluster(*cluster, NULL, 0);
1658 }
1659
1660 FillRow();
1661 fRowCl->GetArray()->Clear("c");
1662
1663 } // for (fRow = 0; fRow < nRows; fRow++) {
1664 if (nClusterSector!=clusterArray->GetEntriesFast()) {
1665 AliError(Form("Failed to read %d out of %d HLT clusters",
1666 clusterArray->GetEntriesFast()-nClusterSector,
1667 clusterArray->GetEntriesFast()));
1668 }
1669 fNclusters+=nClusterSector;
1670 } // for(fSector = 0; fSector < kNS; fSector++) {
03fe1804 1671
6e1c6e9f 1672 pClusterAccess->Clear("event");
03fe1804 1673
53a11d77 1674 Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters);
1675
1676 return 0;
1677}