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