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