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