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