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