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