]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCclustererMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
1. Coding violatation partial fixies
[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");
022ee144 635 }
239f39b1 636 transform->SetCurrentRecoParam((AliTPCRecoParam*)fRecoParam);
022ee144 637 Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
638 Int_t i[1]={fSector};
639 transform->Transform(x,i,0,1);
640 c.SetX(x[0]);
641 c.SetY(x[1]);
642 c.SetZ(x[2]);
643 //
644 //
1c53abe2 645 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
1c53abe2 646 c.SetType(-(c.GetType()+3)); //edge clusters
647 }
648 if (fLoop==2) c.SetType(100);
eea36fac 649 if (!AcceptCluster(&c)) return;
1c53abe2 650
98ee6d31 651 // select output
652 TClonesArray * arr = 0;
653 AliTPCclusterMI * cl = 0;
654
655 if(fBClonesArray==kFALSE) {
656 arr = fRowCl->GetArray();
657 cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
658 } else {
659 cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c);
660 }
661
d101caf3 662 // if (fRecoParam->DumpSignal() &&matrix ) {
663// Int_t nbins=0;
664// Float_t *graph =0;
665// if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
666// nbins = fMaxTime;
667// graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
668// }
669// AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
670// cl->SetInfo(info);
671// }
b127a65f 672 if (!fRecoParam->DumpSignal()) {
673 cl->SetInfo(0);
674 }
eea36fac 675
676 if (AliTPCReconstructor::StreamLevel()>1) {
b7a563f9 677 Float_t xyz[3];
678 cl->GetGlobalXYZ(xyz);
eea36fac 679 (*fDebugStreamer)<<"Clusters"<<
680 "Cl.="<<cl<<
b7a563f9 681 "gx="<<xyz[0]<<
682 "gy="<<xyz[1]<<
683 "gz="<<xyz[2]<<
eea36fac 684 "\n";
685 }
1c53abe2 686
687 fNcluster++;
688}
689
690
691//_____________________________________________________________________________
f8aae377 692void AliTPCclustererMI::Digits2Clusters()
1c53abe2 693{
694 //-----------------------------------------------------------------
695 // This is a simple cluster finder.
696 //-----------------------------------------------------------------
1c53abe2 697
f8aae377 698 if (!fInput) {
699 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 700 return;
701 }
b7a563f9 702 fRecoParam = AliTPCReconstructor::GetRecoParam();
703 if (!fRecoParam){
704 AliFatal("Can not get the reconstruction parameters");
705 }
706 if(AliTPCReconstructor::StreamLevel()>5) {
707 AliInfo("Parameter Dumps");
708 fParam->Dump();
709 fRecoParam->Dump();
710 }
711
13116aec 712 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 713 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1c53abe2 714 AliSimDigits digarr, *dummy=&digarr;
715 fRowDig = dummy;
716 fInput->GetBranch("Segment")->SetAddress(&dummy);
717 Stat_t nentries = fInput->GetEntries();
718
daac2a58 719 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 720
1c53abe2 721 Int_t nclusters = 0;
13116aec 722
1c53abe2 723 for (Int_t n=0; n<nentries; n++) {
724 fInput->GetEvent(n);
508541c7 725 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 726 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
727 continue;
728 }
508541c7 729 Int_t row = fRow;
13116aec 730 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 731 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 732 //
1c53abe2 733
ebe95b38 734 fRowCl->SetID(digarr.GetID());
735 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
f8aae377 736 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 737
738
f8aae377 739 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 740 fZWidth = fParam->GetZWidth();
741 if (fSector < kNIS) {
f8aae377 742 fMaxPad = fParam->GetNPadsLow(row);
1c53abe2 743 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 744 fPadLength = fParam->GetPadPitchLength(fSector,row);
745 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 746 } else {
f8aae377 747 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 748 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 749 fPadLength = fParam->GetPadPitchLength(fSector,row);
750 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 751 }
752
753
754 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 755 fBins =new Float_t[fMaxBin];
5b33a7f5 756 fSigBins =new Int_t[fMaxBin];
757 fNSigBins = 0;
13116aec 758 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 759
760 if (digarr.First()) //MI change
761 do {
13116aec 762 Float_t dig=digarr.CurrentDigit();
f8aae377 763 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 764 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 765 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 766 Int_t bin = i*fMaxTime+j;
d6021085 767 if (gain>0){
768 fBins[bin]=dig/gain;
769 }else{
770 fBins[bin]=0;
771 }
5b33a7f5 772 fSigBins[fNSigBins++]=bin;
1c53abe2 773 } while (digarr.Next());
774 digarr.ExpandTrackBuffer();
775
940ed1f0 776 FindClusters(noiseROC);
ebe95b38 777 FillRow();
e7de6656 778 fRowCl->GetArray()->Clear("C");
88cb7938 779 nclusters+=fNcluster;
98ee6d31 780
5b33a7f5 781 delete[] fBins;
782 delete[] fSigBins;
88cb7938 783 }
b7a563f9 784
19dd5b2f 785 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
f8aae377 786}
787
e7034a8b 788void AliTPCclustererMI::ProcessSectorData(){
7c9528d8 789 //
790 // Process the data for the current sector
791 //
792
793 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
794 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
795 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
796 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
797 //check the presence of the calibration
798 if (!noiseROC ||!pedestalROC ) {
799 AliError(Form("Missing calibration per sector\t%d\n",fSector));
800 return;
801 }
802 Int_t nRows=fParam->GetNRow(fSector);
803 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
804 Int_t zeroSup = fParam->GetZeroSup();
805 // if (calcPedestal) {
806 if (kFALSE ) {
807 for (Int_t iRow = 0; iRow < nRows; iRow++) {
808 Int_t maxPad = fParam->GetNPads(fSector, iRow);
809
810 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
811 //
812 // Temporary fix for data production - !!!! MARIAN
813 // The noise calibration should take mean and RMS - currently the Gaussian fit used
814 // In case of double peak - the pad should be rejected
815 //
816 // Line mean - if more than given digits over threshold - make a noise calculation
817 // and pedestal substration
e7034a8b 818 if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue;
7c9528d8 819 //
e7034a8b 820 if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
821 Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3];
7c9528d8 822 //Float_t pedestal = TMath::Median(fMaxTime, p);
823 Int_t id[3] = {fSector, iRow, iPad-3};
824 // calib values
825 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
826 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
827 Double_t rmsEvent = rmsCalib;
828 Double_t pedestalEvent = pedestalCalib;
829 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
830 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
831 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
832
833 //
834 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
835 Int_t bin = iPad*fMaxTime+iTimeBin;
e7034a8b 836 fAllBins[iRow][bin] -= pedestalEvent;
7c9528d8 837 if (iTimeBin < fRecoParam->GetFirstBin())
e7034a8b 838 fAllBins[iRow][bin] = 0;
7c9528d8 839 if (iTimeBin > fRecoParam->GetLastBin())
e7034a8b 840 fAllBins[iRow][bin] = 0;
841 if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
842 fAllBins[iRow][bin] = 0;
843 if (fAllBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
844 fAllBins[iRow][bin] = 0;
845 if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
7c9528d8 846 }
847 }
848 }
849 }
850
6fbe1e5c 851 if (AliTPCReconstructor::StreamLevel()>5) {
7c9528d8 852 for (Int_t iRow = 0; iRow < nRows; iRow++) {
853 Int_t maxPad = fParam->GetNPads(fSector,iRow);
854
855 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
856 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
857 Int_t bin = iPad*fMaxTime+iTimeBin;
e7034a8b 858 Float_t signal = fAllBins[iRow][bin];
7c9528d8 859 if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
860 Double_t x[]={iRow,iPad-3,iTimeBin-3};
861 Int_t i[]={fSector};
862 AliTPCTransform trafo;
863 trafo.Transform(x,i,0,1);
864 Double_t gx[3]={x[0],x[1],x[2]};
865 trafo.RotatedGlobal2Global(fSector,gx);
e7034a8b 866 // fAllSigBins[iRow][fAllNSigBins[iRow]++]
867 Int_t rowsigBins = fAllNSigBins[iRow];
868 Int_t first=fAllSigBins[iRow][0];
7c9528d8 869 Int_t last= 0;
e7034a8b 870 // if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1];
7c9528d8 871
6fbe1e5c 872 if (AliTPCReconstructor::StreamLevel()>5) {
7c9528d8 873 (*fDebugStreamer)<<"Digits"<<
874 "sec="<<fSector<<
875 "row="<<iRow<<
876 "pad="<<iPad<<
877 "time="<<iTimeBin<<
878 "sig="<<signal<<
879 "x="<<x[0]<<
880 "y="<<x[1]<<
881 "z="<<x[2]<<
882 "gx="<<gx[0]<<
883 "gy="<<gx[1]<<
884 "gz="<<gx[2]<<
885 //
886 "rowsigBins="<<rowsigBins<<
887 "first="<<first<<
888 "last="<<last<<
889 "\n";
890 }
891 }
892 }
893 }
894 }
895 }
896
897 // Now loop over rows and find clusters
898 for (fRow = 0; fRow < nRows; fRow++) {
899 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
900 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
901
902 fRx = fParam->GetPadRowRadii(fSector, fRow);
903 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
904 fPadWidth = fParam->GetPadPitchWidth();
905 fMaxPad = fParam->GetNPads(fSector,fRow);
906 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
907
e7034a8b 908 fBins = fAllBins[fRow];
909 fSigBins = fAllSigBins[fRow];
910 fNSigBins = fAllNSigBins[fRow];
7c9528d8 911
912 FindClusters(noiseROC);
913
914 FillRow();
e7de6656 915 if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C");
7c9528d8 916 fNclusters += fNcluster;
917
918 } // End of loop to find clusters
919}
920
921
f8aae377 922void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
923{
924//-----------------------------------------------------------------
22c352f8 925// This is a cluster finder for the TPC raw data.
926// The method assumes NO ordering of the altro channels.
927// The pedestal subtraction can be switched on and off
928// using an option of the TPC reconstructor
f8aae377 929//-----------------------------------------------------------------
b7a563f9 930 fRecoParam = AliTPCReconstructor::GetRecoParam();
931 if (!fRecoParam){
932 AliFatal("Can not get the reconstruction parameters");
933 }
934 if(AliTPCReconstructor::StreamLevel()>5) {
935 AliInfo("Parameter Dumps");
936 fParam->Dump();
937 fRecoParam->Dump();
938 }
f8aae377 939 fRowDig = NULL;
e7034a8b 940
22c352f8 941 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
7c9528d8 942 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
943 //
944 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
945 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
946 if (fEventHeader){
947 fTimeStamp = fEventHeader->Get("Timestamp");
948 fEventType = fEventHeader->Get("Type");
fefec90f 949 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
950 transform->SetCurrentTimeStamp(fTimeStamp);
951 transform->SetCurrentRun(rawReader->GetRunNumber());
7c9528d8 952 }
953
954 // creaate one TClonesArray for all clusters
955 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
956 // reset counter
957 fNclusters = 0;
958
959 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
960// const Int_t kNIS = fParam->GetNInnerSector();
961// const Int_t kNOS = fParam->GetNOuterSector();
962// const Int_t kNS = kNIS + kNOS;
963 fZWidth = fParam->GetZWidth();
964 Int_t zeroSup = fParam->GetZeroSup();
965 //
e7034a8b 966 // Clean-up
7c9528d8 967 //
e7034a8b 968 AliTPCROC * roc = AliTPCROC::Instance();
7c9528d8 969 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
970 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
7c9528d8 971 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
972 //
973 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
e7034a8b 974 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
975 fAllNSigBins[iRow]=0;
7c9528d8 976 }
977
978 Int_t prevSector=-1;
979 rawReader->Reset();
980 Int_t digCounter=0;
981 //
982 // Loop over DDLs
983 //
984 const Int_t kNIS = fParam->GetNInnerSector();
985 const Int_t kNOS = fParam->GetNOuterSector();
986 const Int_t kNS = kNIS + kNOS;
987
988 for(fSector = 0; fSector < kNS; fSector++) {
989
990 Int_t nRows = 0;
991 Int_t nDDLs = 0, indexDDL = 0;
992 if (fSector < kNIS) {
993 nRows = fParam->GetNRowLow();
994 fSign = (fSector < kNIS/2) ? 1 : -1;
995 nDDLs = 2;
996 indexDDL = fSector * 2;
997 }
998 else {
999 nRows = fParam->GetNRowUp();
1000 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1001 nDDLs = 4;
1002 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1003 }
1004
1005 // load the raw data for corresponding DDLs
1006 rawReader->Reset();
1007 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1008
1009 while (input.NextDDL()){
1010 if (input.GetSector() != fSector)
1011 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1012
fefec90f 1013 //Int_t nRows = fParam->GetNRow(fSector);
7c9528d8 1014
1015 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
1016 // Begin loop over altro data
1017 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1018 Float_t gain =1;
1019
1020 //loop over pads
1021 while ( input.NextChannel() ) {
1022 Int_t iRow = input.GetRow();
1023 if (iRow < 0){
1024 continue;
1025 }
1026 if (iRow >= nRows){
1027 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1028 iRow, 0, nRows -1));
1029 continue;
1030 }
1031 //pad
1032 Int_t iPad = input.GetPad();
1033 if (iPad < 0 || iPad >= nPadsMax) {
1034 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1035 iPad, 0, nPadsMax-1));
1036 continue;
1037 }
1038 gain = gainROC->GetValue(iRow,iPad);
1039 iPad+=3;
1040
1041 //loop over bunches
1042 while ( input.NextBunch() ){
1043 Int_t startTbin = (Int_t)input.GetStartTimeBin();
1044 Int_t bunchlength = (Int_t)input.GetBunchLength();
1045 const UShort_t *sig = input.GetSignals();
1046 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1047 Int_t iTimeBin=startTbin-iTime;
1048 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1049 continue;
1050 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1051 iTimeBin, 0, iTimeBin -1));
1052 }
1053 iTimeBin+=3;
1054 //signal
1055 Float_t signal=(Float_t)sig[iTime];
1056 if (!calcPedestal && signal <= zeroSup) continue;
1057
1058 if (!calcPedestal) {
1059 Int_t bin = iPad*fMaxTime+iTimeBin;
1060 if (gain>0){
e7034a8b 1061 fAllBins[iRow][bin] = signal/gain;
7c9528d8 1062 }else{
e7034a8b 1063 fAllBins[iRow][bin] =0;
7c9528d8 1064 }
e7034a8b 1065 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
7c9528d8 1066 }else{
e7034a8b 1067 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
7c9528d8 1068 }
e7034a8b 1069 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
7c9528d8 1070
1071 // Temporary
1072 digCounter++;
1073 }// end loop signals in bunch
1074 }// end loop bunches
1075 } // end loop pads
1076 //
1077 //
1078 //
1079 //
1080 // Now loop over rows and perform pedestal subtraction
1081 if (digCounter==0) continue;
1082 } // End of loop over sectors
1083 //process last sector
1084 if ( digCounter>0 ){
e7034a8b 1085 ProcessSectorData();
7c9528d8 1086 for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1087 Int_t maxPad = fParam->GetNPads(fSector,iRow);
1088 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
e7034a8b 1089 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1090 fAllNSigBins[iRow] = 0;
7c9528d8 1091 }
1092 prevSector=fSector;
1093 digCounter=0;
1094 }
1095 }
1096
7c9528d8 1097 if (rawReader->GetEventId() && fOutput ){
3a4b0aed 1098 Info("Digits2Clusters", "File %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
7c9528d8 1099 }
1100
1101 if(rawReader->GetEventId()) {
3a4b0aed 1102 Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
7c9528d8 1103 }
1104
1105 if(fBClonesArray) {
1106 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1107 }
1108}
1109
1110
1111
1112
1113
1114void AliTPCclustererMI::Digits2ClustersOld
1115(AliRawReader* rawReader)
1116{
1117//-----------------------------------------------------------------
1118// This is a cluster finder for the TPC raw data.
1119// The method assumes NO ordering of the altro channels.
1120// The pedestal subtraction can be switched on and off
1121// using an option of the TPC reconstructor
1122//-----------------------------------------------------------------
1123 fRecoParam = AliTPCReconstructor::GetRecoParam();
1124 if (!fRecoParam){
1125 AliFatal("Can not get the reconstruction parameters");
1126 }
1127 if(AliTPCReconstructor::StreamLevel()>5) {
1128 AliInfo("Parameter Dumps");
1129 fParam->Dump();
1130 fRecoParam->Dump();
1131 }
1132 fRowDig = NULL;
e7034a8b 1133
7c9528d8 1134 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
d6834f5f 1135 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1136 //
1137 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 1138 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1139 if (fEventHeader){
1140 fTimeStamp = fEventHeader->Get("Timestamp");
1141 fEventType = fEventHeader->Get("Type");
1142 }
1143
98ee6d31 1144 // creaate one TClonesArray for all clusters
44adbd4b 1145 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
98ee6d31 1146 // reset counter
1147 fNclusters = 0;
88cb7938 1148
daac2a58 1149 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 1150 const Int_t kNIS = fParam->GetNInnerSector();
1151 const Int_t kNOS = fParam->GetNOuterSector();
1152 const Int_t kNS = kNIS + kNOS;
1153 fZWidth = fParam->GetZWidth();
1154 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 1155 //
e7034a8b 1156 // Clean-up
adefcaa6 1157 //
e7034a8b 1158
1159 AliTPCROC * roc = AliTPCROC::Instance();
adefcaa6 1160 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1161 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
adefcaa6 1162 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1163 //
1164 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
e7034a8b 1165 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1166 fAllNSigBins[iRow]=0;
adefcaa6 1167 }
1168 //
22c352f8 1169 // Loop over sectors
adefcaa6 1170 //
22c352f8 1171 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 1172
22c352f8 1173 Int_t nRows = 0;
1174 Int_t nDDLs = 0, indexDDL = 0;
1175 if (fSector < kNIS) {
1176 nRows = fParam->GetNRowLow();
1177 fSign = (fSector < kNIS/2) ? 1 : -1;
1178 nDDLs = 2;
1179 indexDDL = fSector * 2;
1180 }
1181 else {
1182 nRows = fParam->GetNRowUp();
1183 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1184 nDDLs = 4;
1185 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1186 }
1187
98ee6d31 1188 // load the raw data for corresponding DDLs
1189 rawReader->Reset();
1190 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1191
1192 // select only good sector
1193 input.Next();
1194 if(input.GetSector() != fSector) continue;
1195
1196 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
98ee6d31 1197
22c352f8 1198 for (Int_t iRow = 0; iRow < nRows; iRow++) {
1199 Int_t maxPad;
1200 if (fSector < kNIS)
7c9528d8 1201 maxPad = fParam->GetNPadsLow(iRow);
22c352f8 1202 else
7c9528d8 1203 maxPad = fParam->GetNPadsUp(iRow);
22c352f8 1204
1205 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
e7034a8b 1206 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1207 fAllNSigBins[iRow] = 0;
22c352f8 1208 }
1209
adefcaa6 1210 Int_t digCounter=0;
22c352f8 1211 // Begin loop over altro data
adefcaa6 1212 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1213 Float_t gain =1;
1214 Int_t lastPad=-1;
98ee6d31 1215
1216 input.Reset();
22c352f8 1217 while (input.Next()) {
22c352f8 1218 if (input.GetSector() != fSector)
7c9528d8 1219 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1220
22c352f8 1221
22c352f8 1222 Int_t iRow = input.GetRow();
9546a7ff 1223 if (iRow < 0){
7c9528d8 1224 continue;
9546a7ff 1225 }
1226
3f0af4ba 1227 if (iRow < 0 || iRow >= nRows){
7c9528d8 1228 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 1229 iRow, 0, nRows -1));
7c9528d8 1230 continue;
3f0af4ba 1231 }
adefcaa6 1232 //pad
1233 Int_t iPad = input.GetPad();
3f0af4ba 1234 if (iPad < 0 || iPad >= nPadsMax) {
7c9528d8 1235 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 1236 iPad, 0, nPadsMax-1));
7c9528d8 1237 continue;
3f0af4ba 1238 }
adefcaa6 1239 if (iPad!=lastPad){
7c9528d8 1240 gain = gainROC->GetValue(iRow,iPad);
1241 lastPad = iPad;
adefcaa6 1242 }
1243 iPad+=3;
1244 //time
1245 Int_t iTimeBin = input.GetTime();
daac2a58 1246 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
7c9528d8 1247 continue;
1248 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 1249 iTimeBin, 0, iTimeBin -1));
daac2a58 1250 }
adefcaa6 1251 iTimeBin+=3;
3f0af4ba 1252
adefcaa6 1253 //signal
22c352f8 1254 Float_t signal = input.GetSignal();
7c9528d8 1255 if (!calcPedestal && signal <= zeroSup) continue;
9546a7ff 1256
940ed1f0 1257 if (!calcPedestal) {
7c9528d8 1258 Int_t bin = iPad*fMaxTime+iTimeBin;
1259 if (gain>0){
e7034a8b 1260 fAllBins[iRow][bin] = signal/gain;
7c9528d8 1261 }else{
e7034a8b 1262 fAllBins[iRow][bin] =0;
7c9528d8 1263 }
e7034a8b 1264 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
940ed1f0 1265 }else{
e7034a8b 1266 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
940ed1f0 1267 }
e7034a8b 1268 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 1269
1270 // Temporary
1271 digCounter++;
22c352f8 1272 } // End of the loop over altro data
adefcaa6 1273 //
1274 //
aba7fdfc 1275 //
1276 //
22c352f8 1277 // Now loop over rows and perform pedestal subtraction
adefcaa6 1278 if (digCounter==0) continue;
e7034a8b 1279 ProcessSectorData();
22c352f8 1280 } // End of loop over sectors
98ee6d31 1281
9546a7ff 1282 if (rawReader->GetEventId() && fOutput ){
98ee6d31 1283 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
9546a7ff 1284 }
ebe95b38 1285
db3576a7 1286 if(rawReader->GetEventId()) {
98ee6d31 1287 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1288 }
1289
1290 if(fBClonesArray) {
1291 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
db3576a7 1292 }
f8aae377 1293}
1294
940ed1f0 1295void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 1296{
940ed1f0 1297
1298 //
1299 // add virtual charge at the edge
1300 //
7041b196 1301 Double_t kMaxDumpSize = 500000;
ebe95b38 1302 if (!fOutput) {
1303 fBDumpSignal =kFALSE;
1304 }else{
1305 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1306 }
1307
f8aae377 1308 fNcluster=0;
f8aae377 1309 fLoop=1;
a1ec4d07 1310 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 1311 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1312 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1313 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1314 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1315 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1316 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
f47588e0 1317 Int_t useOnePadCluster = fRecoParam->GetUseOnePadCluster();
5b33a7f5 1318 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1319 Int_t i = fSigBins[iSig];
1bfaa9e9 1320 if (i%fMaxTime<=crtime) continue;
5b33a7f5 1321 Float_t *b = &fBins[i];
940ed1f0 1322 //absolute custs
03fe1804 1323 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1324 //
f47588e0 1325 if (useOnePadCluster==0){
1326 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1327 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1328 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1329 }
03fe1804 1330 //
1331 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1332 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1333 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1334 //
1335 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1336 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1337 // sigma cuts
1338 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1339 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1340 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1341
d101caf3 1342 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1343 Int_t dummy=0;
1344 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1345
f8aae377 1346 //}
1347 }
8569a2b0 1348}
65c045f0 1349
eea36fac 1350Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1351 //
679026e6 1352 // Currently hack to filter digital noise (15.06.2008)
eea36fac 1353 // To be parameterized in the AliTPCrecoParam
1354 // More inteligent way to be used in future
679026e6 1355 // Acces to the proper pedestal file needed
eea36fac 1356 //
1357 if (cl->GetMax()<400) return kTRUE;
1358 Double_t ratio = cl->GetQ()/cl->GetMax();
679026e6 1359 if (cl->GetMax()>700){
1360 if ((ratio - int(ratio)>0.8)) return kFALSE;
1361 }
eea36fac 1362 if ((ratio - int(ratio)<0.95)) return kTRUE;
1363 return kFALSE;
eea36fac 1364}
1365
65c045f0 1366
940ed1f0 1367Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1368 //
1369 // process signal on given pad - + streaming of additional information in special mode
1370 //
1371 // id[0] - sector
1372 // id[1] - row
adefcaa6 1373 // id[2] - pad
1374
1375 //
1376 // ESTIMATE pedestal and the noise
1377 //
adefcaa6 1378 const Int_t kPedMax = 100;
1379 Float_t max = 0;
1380 Float_t maxPos = 0;
1381 Int_t median = -1;
1382 Int_t count0 = 0;
1383 Int_t count1 = 0;
940ed1f0 1384 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1385 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
b7a563f9 1386 Int_t firstBin = fRecoParam->GetFirstBin();
65c045f0 1387 //
adefcaa6 1388 UShort_t histo[kPedMax];
f174362f 1389 //memset(histo,0,kPedMax*sizeof(UShort_t));
b842afde 1390 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
adefcaa6 1391 for (Int_t i=0; i<fMaxTime; i++){
1392 if (signal[i]<=0) continue;
940ed1f0 1393 if (signal[i]>max && i>firstBin) {
65c045f0 1394 max = signal[i];
1395 maxPos = i;
adefcaa6 1396 }
1397 if (signal[i]>kPedMax-1) continue;
1398 histo[int(signal[i]+0.5)]++;
1399 count0++;
65c045f0 1400 }
7fef31a6 1401 //
adefcaa6 1402 for (Int_t i=1; i<kPedMax; i++){
1403 if (count1<count0*0.5) median=i;
1404 count1+=histo[i];
1405 }
1406 // truncated mean
65c045f0 1407 //
7041b196 1408 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1409 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1410 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1411 //
1412 for (Int_t idelta=1; idelta<10; idelta++){
1413 if (median-idelta<=0) continue;
1414 if (median+idelta>kPedMax) continue;
1415 if (count06<0.6*count1){
1416 count06+=histo[median-idelta];
1417 mean06 +=histo[median-idelta]*(median-idelta);
1418 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1419 count06+=histo[median+idelta];
1420 mean06 +=histo[median+idelta]*(median+idelta);
1421 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1422 }
1423 if (count09<0.9*count1){
1424 count09+=histo[median-idelta];
1425 mean09 +=histo[median-idelta]*(median-idelta);
1426 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1427 count09+=histo[median+idelta];
1428 mean09 +=histo[median+idelta]*(median+idelta);
1429 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1430 }
adefcaa6 1431 if (count10<0.95*count1){
1432 count10+=histo[median-idelta];
1433 mean +=histo[median-idelta]*(median-idelta);
1434 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1435 count10+=histo[median+idelta];
1436 mean +=histo[median+idelta]*(median+idelta);
1437 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1438 }
1439 }
3f0af4ba 1440 if (count10) {
1441 mean /=count10;
1442 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1443 }
1444 if (count06) {
1445 mean06/=count06;
1446 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1447 }
1448 if (count09) {
1449 mean09/=count09;
1450 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1451 }
940ed1f0 1452 rmsEvent = rms09;
adefcaa6 1453 //
940ed1f0 1454 pedestalEvent = median;
adefcaa6 1455 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1456 //
1457 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1458 //
1459 // Dump mean signal info
1460 //
b2426624 1461 if (AliTPCReconstructor::StreamLevel()>0) {
1462 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1463 "TimeStamp="<<fTimeStamp<<
1464 "EventType="<<fEventType<<
adefcaa6 1465 "Sector="<<uid[0]<<
1466 "Row="<<uid[1]<<
1467 "Pad="<<uid[2]<<
1468 "Max="<<max<<
1469 "MaxPos="<<maxPos<<
1470 //
1471 "Median="<<median<<
1472 "Mean="<<mean<<
1473 "RMS="<<rms<<
1474 "Mean06="<<mean06<<
1475 "RMS06="<<rms06<<
1476 "Mean09="<<mean09<<
1477 "RMS09="<<rms09<<
940ed1f0 1478 "RMSCalib="<<rmsCalib<<
1479 "PedCalib="<<pedestalCalib<<
adefcaa6 1480 "\n";
b2426624 1481 }
adefcaa6 1482 //
1483 // fill pedestal histogram
1484 //
65c045f0 1485 //
adefcaa6 1486 //
1487 //
1488 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1489 Float_t *dsignal = new Float_t[nchannels];
1490 Float_t *dtime = new Float_t[nchannels];
1491 for (Int_t i=0; i<nchannels; i++){
1492 dtime[i] = i;
1493 dsignal[i] = signal[i];
1494 }
03fe1804 1495
11441748 1496 TGraph * graph=0;
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]<<
1508 "Graph="<<graph<<
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";
1520 delete graph;
b2426624 1521 }
7fef31a6 1522
65c045f0 1523 delete [] dsignal;
1524 delete [] dtime;
940ed1f0 1525 if (rms06>fRecoParam->GetMaxNoise()) {
1526 pedestalEvent+=1024.;
1527 return 1024+median; // sign noisy channel in debug mode
1528 }
65c045f0 1529 return median;
1530}
1531
1532
1533
03fe1804 1534
03fe1804 1535