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