]> 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
Removing the AliESDCosmicTrack.
[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 }
673
53a11d77 674 //-----------------------------------------------------------------
675 // Use HLT clusters
676 //-----------------------------------------------------------------
49332743 677 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 678 AliInfo("Using HLT clusters for TPC off-line reconstruction");
679 fZWidth = fParam->GetZWidth();
9abc660b 680 Int_t iResult = ReadHLTClusters();
53a11d77 681
9abc660b 682 // HLT clusters present
bbdbe223 683 if (iResult >= 0 && fNclusters > 0)
684 return;
685
9abc660b 686 // HLT clusters not present
bbdbe223 687 if (iResult < 0 || fNclusters == 0) {
688 if (fUseHLTClusters == 3) {
689 AliError("No HLT clusters present, but requested.");
9abc660b 690 return;
691 }
bbdbe223 692 else {
1598ba75 693 AliInfo("Now trying to read from TPC RAW");
bbdbe223 694 }
9abc660b 695 }
696 // Some other problem during cluster reading
697 else {
698 AliWarning("Some problem while unpacking of HLT clusters.");
699 return;
700 }
49332743 701 } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 702
703 //-----------------------------------------------------------------
704 // Run TPC off-line clusterer
705 //-----------------------------------------------------------------
13116aec 706 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 707 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1c53abe2 708 AliSimDigits digarr, *dummy=&digarr;
709 fRowDig = dummy;
710 fInput->GetBranch("Segment")->SetAddress(&dummy);
711 Stat_t nentries = fInput->GetEntries();
712
daac2a58 713 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 714
1c53abe2 715 Int_t nclusters = 0;
13116aec 716
1c53abe2 717 for (Int_t n=0; n<nentries; n++) {
718 fInput->GetEvent(n);
508541c7 719 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 720 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
721 continue;
722 }
508541c7 723 Int_t row = fRow;
13116aec 724 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 725 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 726 //
1c53abe2 727
ebe95b38 728 fRowCl->SetID(digarr.GetID());
729 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
f8aae377 730 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 731
732
f8aae377 733 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 734 fZWidth = fParam->GetZWidth();
735 if (fSector < kNIS) {
f8aae377 736 fMaxPad = fParam->GetNPadsLow(row);
bbdbe223 737 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 738 fPadLength = fParam->GetPadPitchLength(fSector,row);
739 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 740 } else {
f8aae377 741 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 742 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 743 fPadLength = fParam->GetPadPitchLength(fSector,row);
744 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 745 }
746
747
748 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 749 fBins =new Float_t[fMaxBin];
5b33a7f5 750 fSigBins =new Int_t[fMaxBin];
751 fNSigBins = 0;
13116aec 752 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 753
754 if (digarr.First()) //MI change
755 do {
13116aec 756 Float_t dig=digarr.CurrentDigit();
f8aae377 757 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 758 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 759 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 760 Int_t bin = i*fMaxTime+j;
d6021085 761 if (gain>0){
762 fBins[bin]=dig/gain;
763 }else{
764 fBins[bin]=0;
765 }
5b33a7f5 766 fSigBins[fNSigBins++]=bin;
1c53abe2 767 } while (digarr.Next());
768 digarr.ExpandTrackBuffer();
769
940ed1f0 770 FindClusters(noiseROC);
ebe95b38 771 FillRow();
e7de6656 772 fRowCl->GetArray()->Clear("C");
88cb7938 773 nclusters+=fNcluster;
98ee6d31 774
5b33a7f5 775 delete[] fBins;
776 delete[] fSigBins;
88cb7938 777 }
b7a563f9 778
19dd5b2f 779 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
bbdbe223 780
781 if (fUseHLTClusters == 2 && nclusters == 0) {
782 AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
783
784 fZWidth = fParam->GetZWidth();
785 ReadHLTClusters();
786 }
f8aae377 787}
788
e7034a8b 789void AliTPCclustererMI::ProcessSectorData(){
7c9528d8 790 //
791 // Process the data for the current sector
792 //
793
794 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
795 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
796 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
797 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
798 //check the presence of the calibration
799 if (!noiseROC ||!pedestalROC ) {
800 AliError(Form("Missing calibration per sector\t%d\n",fSector));
801 return;
802 }
803 Int_t nRows=fParam->GetNRow(fSector);
804 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
805 Int_t zeroSup = fParam->GetZeroSup();
806 // if (calcPedestal) {
807 if (kFALSE ) {
808 for (Int_t iRow = 0; iRow < nRows; iRow++) {
809 Int_t maxPad = fParam->GetNPads(fSector, iRow);
810
811 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
812 //
813 // Temporary fix for data production - !!!! MARIAN
814 // The noise calibration should take mean and RMS - currently the Gaussian fit used
815 // In case of double peak - the pad should be rejected
816 //
817 // Line mean - if more than given digits over threshold - make a noise calculation
818 // and pedestal substration
e7034a8b 819 if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue;
7c9528d8 820 //
e7034a8b 821 if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
822 Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3];
7c9528d8 823 //Float_t pedestal = TMath::Median(fMaxTime, p);
824 Int_t id[3] = {fSector, iRow, iPad-3};
825 // calib values
826 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
827 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
828 Double_t rmsEvent = rmsCalib;
829 Double_t pedestalEvent = pedestalCalib;
830 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
831 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
832 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
833
834 //
835 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
836 Int_t bin = iPad*fMaxTime+iTimeBin;
e7034a8b 837 fAllBins[iRow][bin] -= pedestalEvent;
7c9528d8 838 if (iTimeBin < fRecoParam->GetFirstBin())
e7034a8b 839 fAllBins[iRow][bin] = 0;
7c9528d8 840 if (iTimeBin > fRecoParam->GetLastBin())
e7034a8b 841 fAllBins[iRow][bin] = 0;
842 if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
843 fAllBins[iRow][bin] = 0;
844 if (fAllBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
845 fAllBins[iRow][bin] = 0;
846 if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
7c9528d8 847 }
848 }
849 }
850 }
851
6fbe1e5c 852 if (AliTPCReconstructor::StreamLevel()>5) {
7c9528d8 853 for (Int_t iRow = 0; iRow < nRows; iRow++) {
854 Int_t maxPad = fParam->GetNPads(fSector,iRow);
855
856 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
857 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
858 Int_t bin = iPad*fMaxTime+iTimeBin;
e7034a8b 859 Float_t signal = fAllBins[iRow][bin];
7c9528d8 860 if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
861 Double_t x[]={iRow,iPad-3,iTimeBin-3};
862 Int_t i[]={fSector};
863 AliTPCTransform trafo;
864 trafo.Transform(x,i,0,1);
865 Double_t gx[3]={x[0],x[1],x[2]};
866 trafo.RotatedGlobal2Global(fSector,gx);
e7034a8b 867 // fAllSigBins[iRow][fAllNSigBins[iRow]++]
868 Int_t rowsigBins = fAllNSigBins[iRow];
869 Int_t first=fAllSigBins[iRow][0];
7c9528d8 870 Int_t last= 0;
e7034a8b 871 // if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1];
7c9528d8 872
6fbe1e5c 873 if (AliTPCReconstructor::StreamLevel()>5) {
7c9528d8 874 (*fDebugStreamer)<<"Digits"<<
875 "sec="<<fSector<<
876 "row="<<iRow<<
877 "pad="<<iPad<<
878 "time="<<iTimeBin<<
879 "sig="<<signal<<
880 "x="<<x[0]<<
881 "y="<<x[1]<<
882 "z="<<x[2]<<
883 "gx="<<gx[0]<<
884 "gy="<<gx[1]<<
885 "gz="<<gx[2]<<
886 //
887 "rowsigBins="<<rowsigBins<<
888 "first="<<first<<
889 "last="<<last<<
890 "\n";
891 }
892 }
893 }
894 }
895 }
896 }
897
898 // Now loop over rows and find clusters
899 for (fRow = 0; fRow < nRows; fRow++) {
900 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
901 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
902
903 fRx = fParam->GetPadRowRadii(fSector, fRow);
904 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
905 fPadWidth = fParam->GetPadPitchWidth();
906 fMaxPad = fParam->GetNPads(fSector,fRow);
907 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
908
e7034a8b 909 fBins = fAllBins[fRow];
910 fSigBins = fAllSigBins[fRow];
911 fNSigBins = fAllNSigBins[fRow];
7c9528d8 912
913 FindClusters(noiseROC);
914
915 FillRow();
e7de6656 916 if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C");
7c9528d8 917 fNclusters += fNcluster;
918
919 } // End of loop to find clusters
920}
921
922
f8aae377 923void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
924{
925//-----------------------------------------------------------------
22c352f8 926// This is a cluster finder for the TPC raw data.
927// The method assumes NO ordering of the altro channels.
928// The pedestal subtraction can be switched on and off
929// using an option of the TPC reconstructor
f8aae377 930//-----------------------------------------------------------------
b7a563f9 931 fRecoParam = AliTPCReconstructor::GetRecoParam();
932 if (!fRecoParam){
933 AliFatal("Can not get the reconstruction parameters");
934 }
935 if(AliTPCReconstructor::StreamLevel()>5) {
936 AliInfo("Parameter Dumps");
937 fParam->Dump();
938 fRecoParam->Dump();
939 }
f8aae377 940 fRowDig = NULL;
e7034a8b 941
53a11d77 942 //-----------------------------------------------------------------
943 // Use HLT clusters
944 //-----------------------------------------------------------------
49332743 945 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 946 AliInfo("Using HLT clusters for TPC off-line reconstruction");
947 fZWidth = fParam->GetZWidth();
9abc660b 948 Int_t iResult = ReadHLTClusters();
53a11d77 949
9abc660b 950 // HLT clusters present
bbdbe223 951 if (iResult >= 0 && fNclusters > 0)
9abc660b 952 return;
bbdbe223 953
9abc660b 954 // HLT clusters not present
bbdbe223 955 if (iResult < 0 || fNclusters == 0) {
956 if (fUseHLTClusters == 3) {
957 AliError("No HLT clusters present, but requested.");
9abc660b 958 return;
959 }
bbdbe223 960 else {
961 AliInfo("Now trying to read TPC RAW");
962 }
9abc660b 963 }
964 // Some other problem during cluster reading
965 else {
966 AliWarning("Some problem while unpacking of HLT clusters.");
967 return;
968 }
49332743 969 } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
53a11d77 970
971 //-----------------------------------------------------------------
972 // Run TPC off-line clusterer
973 //-----------------------------------------------------------------
22c352f8 974 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
7c9528d8 975 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
976 //
977 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
978 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
979 if (fEventHeader){
980 fTimeStamp = fEventHeader->Get("Timestamp");
981 fEventType = fEventHeader->Get("Type");
fefec90f 982 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
983 transform->SetCurrentTimeStamp(fTimeStamp);
984 transform->SetCurrentRun(rawReader->GetRunNumber());
7c9528d8 985 }
986
987 // creaate one TClonesArray for all clusters
988 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
989 // reset counter
990 fNclusters = 0;
991
992 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
993// const Int_t kNIS = fParam->GetNInnerSector();
994// const Int_t kNOS = fParam->GetNOuterSector();
995// const Int_t kNS = kNIS + kNOS;
996 fZWidth = fParam->GetZWidth();
997 Int_t zeroSup = fParam->GetZeroSup();
998 //
e7034a8b 999 // Clean-up
7c9528d8 1000 //
e7034a8b 1001 AliTPCROC * roc = AliTPCROC::Instance();
7c9528d8 1002 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1003 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
7c9528d8 1004 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1005 //
1006 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
e7034a8b 1007 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1008 fAllNSigBins[iRow]=0;
7c9528d8 1009 }
1010
1011 Int_t prevSector=-1;
1012 rawReader->Reset();
1013 Int_t digCounter=0;
1014 //
1015 // Loop over DDLs
1016 //
1017 const Int_t kNIS = fParam->GetNInnerSector();
1018 const Int_t kNOS = fParam->GetNOuterSector();
1019 const Int_t kNS = kNIS + kNOS;
1020
1021 for(fSector = 0; fSector < kNS; fSector++) {
1022
1023 Int_t nRows = 0;
1024 Int_t nDDLs = 0, indexDDL = 0;
1025 if (fSector < kNIS) {
1026 nRows = fParam->GetNRowLow();
1027 fSign = (fSector < kNIS/2) ? 1 : -1;
1028 nDDLs = 2;
1029 indexDDL = fSector * 2;
1030 }
1031 else {
1032 nRows = fParam->GetNRowUp();
1033 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1034 nDDLs = 4;
1035 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1036 }
1037
1038 // load the raw data for corresponding DDLs
1039 rawReader->Reset();
1040 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1041
1042 while (input.NextDDL()){
1043 if (input.GetSector() != fSector)
1044 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1045
fefec90f 1046 //Int_t nRows = fParam->GetNRow(fSector);
7c9528d8 1047
1048 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
1049 // Begin loop over altro data
1050 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1051 Float_t gain =1;
1052
1053 //loop over pads
1054 while ( input.NextChannel() ) {
1055 Int_t iRow = input.GetRow();
1056 if (iRow < 0){
1057 continue;
1058 }
1059 if (iRow >= nRows){
1060 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1061 iRow, 0, nRows -1));
1062 continue;
1063 }
1064 //pad
1065 Int_t iPad = input.GetPad();
1066 if (iPad < 0 || iPad >= nPadsMax) {
1067 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1068 iPad, 0, nPadsMax-1));
1069 continue;
1070 }
1071 gain = gainROC->GetValue(iRow,iPad);
1072 iPad+=3;
1073
1074 //loop over bunches
1075 while ( input.NextBunch() ){
1076 Int_t startTbin = (Int_t)input.GetStartTimeBin();
1077 Int_t bunchlength = (Int_t)input.GetBunchLength();
1078 const UShort_t *sig = input.GetSignals();
1079 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1080 Int_t iTimeBin=startTbin-iTime;
1081 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1082 continue;
1083 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1084 iTimeBin, 0, iTimeBin -1));
1085 }
1086 iTimeBin+=3;
1087 //signal
1088 Float_t signal=(Float_t)sig[iTime];
1089 if (!calcPedestal && signal <= zeroSup) continue;
1090
1091 if (!calcPedestal) {
1092 Int_t bin = iPad*fMaxTime+iTimeBin;
1093 if (gain>0){
e7034a8b 1094 fAllBins[iRow][bin] = signal/gain;
7c9528d8 1095 }else{
e7034a8b 1096 fAllBins[iRow][bin] =0;
7c9528d8 1097 }
e7034a8b 1098 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
7c9528d8 1099 }else{
e7034a8b 1100 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
7c9528d8 1101 }
e7034a8b 1102 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
7c9528d8 1103
1104 // Temporary
1105 digCounter++;
1106 }// end loop signals in bunch
1107 }// end loop bunches
1108 } // end loop pads
1109 //
1110 //
1111 //
1112 //
1113 // Now loop over rows and perform pedestal subtraction
1114 if (digCounter==0) continue;
1115 } // End of loop over sectors
1116 //process last sector
1117 if ( digCounter>0 ){
e7034a8b 1118 ProcessSectorData();
7c9528d8 1119 for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1120 Int_t maxPad = fParam->GetNPads(fSector,iRow);
1121 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
e7034a8b 1122 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1123 fAllNSigBins[iRow] = 0;
7c9528d8 1124 }
1125 prevSector=fSector;
1126 digCounter=0;
1127 }
1128 }
1129
7c9528d8 1130 if (rawReader->GetEventId() && fOutput ){
3a4b0aed 1131 Info("Digits2Clusters", "File %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
7c9528d8 1132 }
1133
1134 if(rawReader->GetEventId()) {
3a4b0aed 1135 Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
7c9528d8 1136 }
1137
1138 if(fBClonesArray) {
1139 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1140 }
bbdbe223 1141
1142 if (fUseHLTClusters == 2 && fNclusters == 0) {
1143 AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
1144
1145 fZWidth = fParam->GetZWidth();
1146 ReadHLTClusters();
1147 }
7c9528d8 1148}
1149
1150
1151
1152
1153
1154void AliTPCclustererMI::Digits2ClustersOld
1155(AliRawReader* rawReader)
1156{
1157//-----------------------------------------------------------------
1158// This is a cluster finder for the TPC raw data.
1159// The method assumes NO ordering of the altro channels.
1160// The pedestal subtraction can be switched on and off
1161// using an option of the TPC reconstructor
1162//-----------------------------------------------------------------
1163 fRecoParam = AliTPCReconstructor::GetRecoParam();
1164 if (!fRecoParam){
1165 AliFatal("Can not get the reconstruction parameters");
1166 }
1167 if(AliTPCReconstructor::StreamLevel()>5) {
1168 AliInfo("Parameter Dumps");
1169 fParam->Dump();
1170 fRecoParam->Dump();
1171 }
1172 fRowDig = NULL;
e7034a8b 1173
7c9528d8 1174 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
d6834f5f 1175 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1176 //
1177 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 1178 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1179 if (fEventHeader){
1180 fTimeStamp = fEventHeader->Get("Timestamp");
1181 fEventType = fEventHeader->Get("Type");
1182 }
1183
98ee6d31 1184 // creaate one TClonesArray for all clusters
44adbd4b 1185 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
98ee6d31 1186 // reset counter
1187 fNclusters = 0;
88cb7938 1188
daac2a58 1189 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 1190 const Int_t kNIS = fParam->GetNInnerSector();
1191 const Int_t kNOS = fParam->GetNOuterSector();
1192 const Int_t kNS = kNIS + kNOS;
1193 fZWidth = fParam->GetZWidth();
1194 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 1195 //
e7034a8b 1196 // Clean-up
adefcaa6 1197 //
e7034a8b 1198
1199 AliTPCROC * roc = AliTPCROC::Instance();
adefcaa6 1200 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1201 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
adefcaa6 1202 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1203 //
1204 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
e7034a8b 1205 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1206 fAllNSigBins[iRow]=0;
adefcaa6 1207 }
1208 //
22c352f8 1209 // Loop over sectors
adefcaa6 1210 //
22c352f8 1211 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 1212
22c352f8 1213 Int_t nRows = 0;
1214 Int_t nDDLs = 0, indexDDL = 0;
1215 if (fSector < kNIS) {
1216 nRows = fParam->GetNRowLow();
1217 fSign = (fSector < kNIS/2) ? 1 : -1;
1218 nDDLs = 2;
1219 indexDDL = fSector * 2;
1220 }
1221 else {
1222 nRows = fParam->GetNRowUp();
1223 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1224 nDDLs = 4;
1225 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1226 }
1227
98ee6d31 1228 // load the raw data for corresponding DDLs
1229 rawReader->Reset();
1230 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1231
1232 // select only good sector
35e039dc 1233 if (!input.Next()) continue;
98ee6d31 1234 if(input.GetSector() != fSector) continue;
1235
1236 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
98ee6d31 1237
22c352f8 1238 for (Int_t iRow = 0; iRow < nRows; iRow++) {
1239 Int_t maxPad;
1240 if (fSector < kNIS)
7c9528d8 1241 maxPad = fParam->GetNPadsLow(iRow);
22c352f8 1242 else
7c9528d8 1243 maxPad = fParam->GetNPadsUp(iRow);
22c352f8 1244
1245 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
e7034a8b 1246 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1247 fAllNSigBins[iRow] = 0;
22c352f8 1248 }
1249
adefcaa6 1250 Int_t digCounter=0;
22c352f8 1251 // Begin loop over altro data
adefcaa6 1252 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1253 Float_t gain =1;
1254 Int_t lastPad=-1;
98ee6d31 1255
1256 input.Reset();
22c352f8 1257 while (input.Next()) {
22c352f8 1258 if (input.GetSector() != fSector)
7c9528d8 1259 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1260
22c352f8 1261
22c352f8 1262 Int_t iRow = input.GetRow();
9546a7ff 1263 if (iRow < 0){
7c9528d8 1264 continue;
9546a7ff 1265 }
1266
3f0af4ba 1267 if (iRow < 0 || iRow >= nRows){
7c9528d8 1268 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 1269 iRow, 0, nRows -1));
7c9528d8 1270 continue;
3f0af4ba 1271 }
adefcaa6 1272 //pad
1273 Int_t iPad = input.GetPad();
3f0af4ba 1274 if (iPad < 0 || iPad >= nPadsMax) {
7c9528d8 1275 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 1276 iPad, 0, nPadsMax-1));
7c9528d8 1277 continue;
3f0af4ba 1278 }
adefcaa6 1279 if (iPad!=lastPad){
7c9528d8 1280 gain = gainROC->GetValue(iRow,iPad);
1281 lastPad = iPad;
adefcaa6 1282 }
1283 iPad+=3;
1284 //time
1285 Int_t iTimeBin = input.GetTime();
daac2a58 1286 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
7c9528d8 1287 continue;
1288 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 1289 iTimeBin, 0, iTimeBin -1));
daac2a58 1290 }
adefcaa6 1291 iTimeBin+=3;
3f0af4ba 1292
adefcaa6 1293 //signal
22c352f8 1294 Float_t signal = input.GetSignal();
7c9528d8 1295 if (!calcPedestal && signal <= zeroSup) continue;
9546a7ff 1296
940ed1f0 1297 if (!calcPedestal) {
7c9528d8 1298 Int_t bin = iPad*fMaxTime+iTimeBin;
1299 if (gain>0){
e7034a8b 1300 fAllBins[iRow][bin] = signal/gain;
7c9528d8 1301 }else{
e7034a8b 1302 fAllBins[iRow][bin] =0;
7c9528d8 1303 }
e7034a8b 1304 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
940ed1f0 1305 }else{
e7034a8b 1306 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
940ed1f0 1307 }
e7034a8b 1308 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 1309
1310 // Temporary
1311 digCounter++;
22c352f8 1312 } // End of the loop over altro data
adefcaa6 1313 //
1314 //
aba7fdfc 1315 //
1316 //
22c352f8 1317 // Now loop over rows and perform pedestal subtraction
adefcaa6 1318 if (digCounter==0) continue;
e7034a8b 1319 ProcessSectorData();
22c352f8 1320 } // End of loop over sectors
98ee6d31 1321
9546a7ff 1322 if (rawReader->GetEventId() && fOutput ){
98ee6d31 1323 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
9546a7ff 1324 }
ebe95b38 1325
db3576a7 1326 if(rawReader->GetEventId()) {
98ee6d31 1327 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1328 }
1329
1330 if(fBClonesArray) {
1331 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
db3576a7 1332 }
f8aae377 1333}
1334
940ed1f0 1335void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 1336{
940ed1f0 1337
1338 //
1339 // add virtual charge at the edge
1340 //
7041b196 1341 Double_t kMaxDumpSize = 500000;
ebe95b38 1342 if (!fOutput) {
1343 fBDumpSignal =kFALSE;
1344 }else{
1345 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1346 }
1347
f8aae377 1348 fNcluster=0;
f8aae377 1349 fLoop=1;
a1ec4d07 1350 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 1351 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1352 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1353 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1354 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1355 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1356 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
f47588e0 1357 Int_t useOnePadCluster = fRecoParam->GetUseOnePadCluster();
5b33a7f5 1358 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1359 Int_t i = fSigBins[iSig];
1bfaa9e9 1360 if (i%fMaxTime<=crtime) continue;
5b33a7f5 1361 Float_t *b = &fBins[i];
940ed1f0 1362 //absolute custs
03fe1804 1363 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1364 //
f47588e0 1365 if (useOnePadCluster==0){
1366 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1367 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1368 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1369 }
03fe1804 1370 //
1371 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1372 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1373 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1374 //
1375 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1376 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1377 // sigma cuts
1378 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1379 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1380 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1381
d101caf3 1382 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1383 Int_t dummy=0;
1384 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1385
f8aae377 1386 //}
1387 }
8569a2b0 1388}
65c045f0 1389
eea36fac 1390Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
80efd770 1391 // -- Depricated --
679026e6 1392 // Currently hack to filter digital noise (15.06.2008)
eea36fac 1393 // To be parameterized in the AliTPCrecoParam
1394 // More inteligent way to be used in future
679026e6 1395 // Acces to the proper pedestal file needed
eea36fac 1396 //
1397 if (cl->GetMax()<400) return kTRUE;
1398 Double_t ratio = cl->GetQ()/cl->GetMax();
679026e6 1399 if (cl->GetMax()>700){
1400 if ((ratio - int(ratio)>0.8)) return kFALSE;
1401 }
eea36fac 1402 if ((ratio - int(ratio)<0.95)) return kTRUE;
1403 return kFALSE;
eea36fac 1404}
1405
65c045f0 1406
940ed1f0 1407Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1408 //
1409 // process signal on given pad - + streaming of additional information in special mode
1410 //
1411 // id[0] - sector
1412 // id[1] - row
adefcaa6 1413 // id[2] - pad
1414
1415 //
1416 // ESTIMATE pedestal and the noise
1417 //
adefcaa6 1418 const Int_t kPedMax = 100;
1419 Float_t max = 0;
1420 Float_t maxPos = 0;
1421 Int_t median = -1;
1422 Int_t count0 = 0;
1423 Int_t count1 = 0;
940ed1f0 1424 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1425 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
b7a563f9 1426 Int_t firstBin = fRecoParam->GetFirstBin();
65c045f0 1427 //
adefcaa6 1428 UShort_t histo[kPedMax];
f174362f 1429 //memset(histo,0,kPedMax*sizeof(UShort_t));
b842afde 1430 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
adefcaa6 1431 for (Int_t i=0; i<fMaxTime; i++){
1432 if (signal[i]<=0) continue;
940ed1f0 1433 if (signal[i]>max && i>firstBin) {
65c045f0 1434 max = signal[i];
1435 maxPos = i;
adefcaa6 1436 }
1437 if (signal[i]>kPedMax-1) continue;
1438 histo[int(signal[i]+0.5)]++;
1439 count0++;
65c045f0 1440 }
7fef31a6 1441 //
adefcaa6 1442 for (Int_t i=1; i<kPedMax; i++){
1443 if (count1<count0*0.5) median=i;
1444 count1+=histo[i];
1445 }
1446 // truncated mean
65c045f0 1447 //
7041b196 1448 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1449 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1450 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1451 //
1452 for (Int_t idelta=1; idelta<10; idelta++){
1453 if (median-idelta<=0) continue;
1454 if (median+idelta>kPedMax) continue;
1455 if (count06<0.6*count1){
1456 count06+=histo[median-idelta];
1457 mean06 +=histo[median-idelta]*(median-idelta);
1458 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1459 count06+=histo[median+idelta];
1460 mean06 +=histo[median+idelta]*(median+idelta);
1461 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1462 }
1463 if (count09<0.9*count1){
1464 count09+=histo[median-idelta];
1465 mean09 +=histo[median-idelta]*(median-idelta);
1466 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1467 count09+=histo[median+idelta];
1468 mean09 +=histo[median+idelta]*(median+idelta);
1469 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1470 }
adefcaa6 1471 if (count10<0.95*count1){
1472 count10+=histo[median-idelta];
1473 mean +=histo[median-idelta]*(median-idelta);
1474 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1475 count10+=histo[median+idelta];
1476 mean +=histo[median+idelta]*(median+idelta);
1477 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1478 }
1479 }
3f0af4ba 1480 if (count10) {
1481 mean /=count10;
1482 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1483 }
1484 if (count06) {
1485 mean06/=count06;
1486 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1487 }
1488 if (count09) {
1489 mean09/=count09;
1490 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1491 }
940ed1f0 1492 rmsEvent = rms09;
adefcaa6 1493 //
940ed1f0 1494 pedestalEvent = median;
adefcaa6 1495 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1496 //
1497 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1498 //
1499 // Dump mean signal info
1500 //
b2426624 1501 if (AliTPCReconstructor::StreamLevel()>0) {
1502 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1503 "TimeStamp="<<fTimeStamp<<
1504 "EventType="<<fEventType<<
adefcaa6 1505 "Sector="<<uid[0]<<
1506 "Row="<<uid[1]<<
1507 "Pad="<<uid[2]<<
1508 "Max="<<max<<
1509 "MaxPos="<<maxPos<<
1510 //
1511 "Median="<<median<<
1512 "Mean="<<mean<<
1513 "RMS="<<rms<<
1514 "Mean06="<<mean06<<
1515 "RMS06="<<rms06<<
1516 "Mean09="<<mean09<<
1517 "RMS09="<<rms09<<
940ed1f0 1518 "RMSCalib="<<rmsCalib<<
1519 "PedCalib="<<pedestalCalib<<
adefcaa6 1520 "\n";
b2426624 1521 }
adefcaa6 1522 //
1523 // fill pedestal histogram
1524 //
65c045f0 1525 //
adefcaa6 1526 //
1527 //
1528 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1529 Float_t *dsignal = new Float_t[nchannels];
1530 Float_t *dtime = new Float_t[nchannels];
1531 for (Int_t i=0; i<nchannels; i++){
1532 dtime[i] = i;
1533 dsignal[i] = signal[i];
1534 }
03fe1804 1535
7fef31a6 1536 //
a525cc34 1537 // Big signals dumping
1538 //
b2426624 1539 if (AliTPCReconstructor::StreamLevel()>0) {
b7a563f9 1540 if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin())
a525cc34 1541 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1542 "TimeStamp="<<fTimeStamp<<
1543 "EventType="<<fEventType<<
1544 "Sector="<<uid[0]<<
1545 "Row="<<uid[1]<<
1546 "Pad="<<uid[2]<<
35e039dc 1547 // "Graph="<<graph<<
a525cc34 1548 "Max="<<max<<
1549 "MaxPos="<<maxPos<<
1550 //
1551 "Median="<<median<<
1552 "Mean="<<mean<<
1553 "RMS="<<rms<<
1554 "Mean06="<<mean06<<
1555 "RMS06="<<rms06<<
1556 "Mean09="<<mean09<<
1557 "RMS09="<<rms09<<
1558 "\n";
b2426624 1559 }
7fef31a6 1560
65c045f0 1561 delete [] dsignal;
1562 delete [] dtime;
940ed1f0 1563 if (rms06>fRecoParam->GetMaxNoise()) {
1564 pedestalEvent+=1024.;
1565 return 1024+median; // sign noisy channel in debug mode
1566 }
65c045f0 1567 return median;
1568}
1569
53a11d77 1570Int_t AliTPCclustererMI::ReadHLTClusters()
1571{
1572 //
1573 // read HLT clusters instead of off line custers,
1574 // used in Digits2Clusters
1575 //
1576
6e1c6e9f 1577 if (!fHLTClusterAccess) {
53a11d77 1578 TClass* pCl=NULL;
1579 ROOT::NewFunc_t pNewFunc=NULL;
1580 do {
1581 pCl=TClass::GetClass("AliHLTTPCClusterAccessHLTOUT");
1582 } while (!pCl && gSystem->Load("libAliHLTTPC.so")==0);
1583 if (!pCl || (pNewFunc=pCl->GetNew())==NULL) {
1584 AliError("can not load class description of AliHLTTPCClusterAccessHLTOUT, aborting ...");
1585 return -1;
1586 }
1587
1588 void* p=(*pNewFunc)(NULL);
1589 if (!p) {
1590 AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT");
1591 return -2;
1592 }
6e1c6e9f 1593 fHLTClusterAccess=reinterpret_cast<TObject*>(p);
53a11d77 1594 }
1595
6e1c6e9f 1596 TObject* pClusterAccess=fHLTClusterAccess;
1597
53a11d77 1598 const Int_t kNIS = fParam->GetNInnerSector();
1599 const Int_t kNOS = fParam->GetNOuterSector();
1600 const Int_t kNS = kNIS + kNOS;
1601 fNclusters = 0;
1602
6e1c6e9f 1603 // make sure that all clusters from the previous event are cleared
1604 pClusterAccess->Clear("event");
53a11d77 1605 for(fSector = 0; fSector < kNS; fSector++) {
1606
bbdbe223 1607 Int_t iResult = 1;
53a11d77 1608 TString param("sector="); param+=fSector;
6e1c6e9f 1609 // prepare for next sector
1610 pClusterAccess->Clear("sector");
bbdbe223 1611 pClusterAccess->Execute("read", param, &iResult);
1612 if (iResult < 0) {
1613 return iResult;
1614 AliError("HLT Clusters can not be found");
1615 }
1616
53a11d77 1617 if (pClusterAccess->FindObject("clusterarray")==NULL) {
1618 AliError("HLT clusters requested, but not cluster array not present");
1619 return -4;
1620 }
1621
1622 TClonesArray* clusterArray=dynamic_cast<TClonesArray*>(pClusterAccess->FindObject("clusterarray"));
1623 if (!clusterArray) {
1624 AliError("HLT cluster array is not of class type TClonesArray");
1625 return -5;
1626 }
1627
1628 AliDebug(4,Form("Reading %d clusters from HLT for sector %d", clusterArray->GetEntriesFast(), fSector));
65c045f0 1629
53a11d77 1630 Int_t nClusterSector=0;
1631 Int_t nRows=fParam->GetNRow(fSector);
65c045f0 1632
53a11d77 1633 for (fRow = 0; fRow < nRows; fRow++) {
1634 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
1635 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
1636 fNcluster=0; // reset clusters per row
1637
1638 fRx = fParam->GetPadRowRadii(fSector, fRow);
1639 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
1640 fPadWidth = fParam->GetPadPitchWidth();
1641 fMaxPad = fParam->GetNPads(fSector,fRow);
1642 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
1643
1644 fBins = fAllBins[fRow];
1645 fSigBins = fAllSigBins[fRow];
1646 fNSigBins = fAllNSigBins[fRow];
1647
1648 for (Int_t i=0; i<clusterArray->GetEntriesFast(); i++) {
1649 if (!clusterArray->At(i))
1650 continue;
1651
1652 AliTPCclusterMI* cluster=dynamic_cast<AliTPCclusterMI*>(clusterArray->At(i));
1653 if (!cluster) continue;
1654 if (cluster->GetRow()!=fRow) continue;
1655 nClusterSector++;
1656 AddCluster(*cluster, NULL, 0);
1657 }
1658
1659 FillRow();
1660 fRowCl->GetArray()->Clear("c");
1661
1662 } // for (fRow = 0; fRow < nRows; fRow++) {
1663 if (nClusterSector!=clusterArray->GetEntriesFast()) {
1664 AliError(Form("Failed to read %d out of %d HLT clusters",
1665 clusterArray->GetEntriesFast()-nClusterSector,
1666 clusterArray->GetEntriesFast()));
1667 }
1668 fNclusters+=nClusterSector;
1669 } // for(fSector = 0; fSector < kNS; fSector++) {
03fe1804 1670
6e1c6e9f 1671 pClusterAccess->Clear("event");
03fe1804 1672
53a11d77 1673 Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters);
1674
1675 return 0;
1676}