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