]> 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
Changes of Jens Wiechula
[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");
913 }
914
915 // creaate one TClonesArray for all clusters
916 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
917 // reset counter
918 fNclusters = 0;
919
920 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
921// const Int_t kNIS = fParam->GetNInnerSector();
922// const Int_t kNOS = fParam->GetNOuterSector();
923// const Int_t kNS = kNIS + kNOS;
924 fZWidth = fParam->GetZWidth();
925 Int_t zeroSup = fParam->GetZeroSup();
926 //
927 //alocate memory for sector - maximal case
928 //
929 Float_t** allBins = NULL;
930 Int_t** allSigBins = NULL;
931 Int_t* allNSigBins = NULL;
932 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
933 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
934 allBins = new Float_t*[nRowsMax];
935 allSigBins = new Int_t*[nRowsMax];
936 allNSigBins = new Int_t[nRowsMax];
937 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
938 //
939 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
940 allBins[iRow] = new Float_t[maxBin];
941 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
942 allSigBins[iRow] = new Int_t[maxBin];
943 allNSigBins[iRow]=0;
944 }
945
946 Int_t prevSector=-1;
947 rawReader->Reset();
948 Int_t digCounter=0;
949 //
950 // Loop over DDLs
951 //
952 const Int_t kNIS = fParam->GetNInnerSector();
953 const Int_t kNOS = fParam->GetNOuterSector();
954 const Int_t kNS = kNIS + kNOS;
955
956 for(fSector = 0; fSector < kNS; fSector++) {
957
958 Int_t nRows = 0;
959 Int_t nDDLs = 0, indexDDL = 0;
960 if (fSector < kNIS) {
961 nRows = fParam->GetNRowLow();
962 fSign = (fSector < kNIS/2) ? 1 : -1;
963 nDDLs = 2;
964 indexDDL = fSector * 2;
965 }
966 else {
967 nRows = fParam->GetNRowUp();
968 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
969 nDDLs = 4;
970 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
971 }
972
973 // load the raw data for corresponding DDLs
974 rawReader->Reset();
975 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
976
977 while (input.NextDDL()){
978 if (input.GetSector() != fSector)
979 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
980
981 Int_t nRows = fParam->GetNRow(fSector);
982
983 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
984 // Begin loop over altro data
985 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
986 Float_t gain =1;
987
988 //loop over pads
989 while ( input.NextChannel() ) {
990 Int_t iRow = input.GetRow();
991 if (iRow < 0){
992 continue;
993 }
994 if (iRow >= nRows){
995 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
996 iRow, 0, nRows -1));
997 continue;
998 }
999 //pad
1000 Int_t iPad = input.GetPad();
1001 if (iPad < 0 || iPad >= nPadsMax) {
1002 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1003 iPad, 0, nPadsMax-1));
1004 continue;
1005 }
1006 gain = gainROC->GetValue(iRow,iPad);
1007 iPad+=3;
1008
1009 //loop over bunches
1010 while ( input.NextBunch() ){
1011 Int_t startTbin = (Int_t)input.GetStartTimeBin();
1012 Int_t bunchlength = (Int_t)input.GetBunchLength();
1013 const UShort_t *sig = input.GetSignals();
1014 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1015 Int_t iTimeBin=startTbin-iTime;
1016 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1017 continue;
1018 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1019 iTimeBin, 0, iTimeBin -1));
1020 }
1021 iTimeBin+=3;
1022 //signal
1023 Float_t signal=(Float_t)sig[iTime];
1024 if (!calcPedestal && signal <= zeroSup) continue;
1025
1026 if (!calcPedestal) {
1027 Int_t bin = iPad*fMaxTime+iTimeBin;
1028 if (gain>0){
1029 allBins[iRow][bin] = signal/gain;
1030 }else{
1031 allBins[iRow][bin] =0;
1032 }
1033 allSigBins[iRow][allNSigBins[iRow]++] = bin;
1034 }else{
1035 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1036 }
1037 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
1038
1039 // Temporary
1040 digCounter++;
1041 }// end loop signals in bunch
1042 }// end loop bunches
1043 } // end loop pads
1044 //
1045 //
1046 //
1047 //
1048 // Now loop over rows and perform pedestal subtraction
1049 if (digCounter==0) continue;
1050 } // End of loop over sectors
1051 //process last sector
1052 if ( digCounter>0 ){
1053 ProcessSectorData(allBins, allSigBins, allNSigBins);
1054 for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1055 Int_t maxPad = fParam->GetNPads(fSector,iRow);
1056 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
1057 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
1058 allNSigBins[iRow] = 0;
1059 }
1060 prevSector=fSector;
1061 digCounter=0;
1062 }
1063 }
1064
1065
1066
1067 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1068 delete [] allBins[iRow];
1069 delete [] allSigBins[iRow];
1070 }
1071 delete [] allBins;
1072 delete [] allSigBins;
1073 delete [] allNSigBins;
1074
1075 if (rawReader->GetEventId() && fOutput ){
1076 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1077 }
1078
1079 if(rawReader->GetEventId()) {
1080 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1081 }
1082
1083 if(fBClonesArray) {
1084 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1085 }
1086}
1087
1088
1089
1090
1091
1092void AliTPCclustererMI::Digits2ClustersOld
1093(AliRawReader* rawReader)
1094{
1095//-----------------------------------------------------------------
1096// This is a cluster finder for the TPC raw data.
1097// The method assumes NO ordering of the altro channels.
1098// The pedestal subtraction can be switched on and off
1099// using an option of the TPC reconstructor
1100//-----------------------------------------------------------------
1101 fRecoParam = AliTPCReconstructor::GetRecoParam();
1102 if (!fRecoParam){
1103 AliFatal("Can not get the reconstruction parameters");
1104 }
1105 if(AliTPCReconstructor::StreamLevel()>5) {
1106 AliInfo("Parameter Dumps");
1107 fParam->Dump();
1108 fRecoParam->Dump();
1109 }
1110 fRowDig = NULL;
1111 AliTPCROC * roc = AliTPCROC::Instance();
1112 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
d6834f5f 1113 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1114 //
1115 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 1116 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1117 if (fEventHeader){
1118 fTimeStamp = fEventHeader->Get("Timestamp");
1119 fEventType = fEventHeader->Get("Type");
1120 }
1121
98ee6d31 1122 // creaate one TClonesArray for all clusters
44adbd4b 1123 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
98ee6d31 1124 // reset counter
1125 fNclusters = 0;
88cb7938 1126
daac2a58 1127 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 1128 const Int_t kNIS = fParam->GetNInnerSector();
1129 const Int_t kNOS = fParam->GetNOuterSector();
1130 const Int_t kNS = kNIS + kNOS;
1131 fZWidth = fParam->GetZWidth();
1132 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 1133 //
1134 //alocate memory for sector - maximal case
1135 //
22c352f8 1136 Float_t** allBins = NULL;
5b33a7f5 1137 Int_t** allSigBins = NULL;
1138 Int_t* allNSigBins = NULL;
adefcaa6 1139 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1140 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1141 allBins = new Float_t*[nRowsMax];
5b33a7f5 1142 allSigBins = new Int_t*[nRowsMax];
1143 allNSigBins = new Int_t[nRowsMax];
adefcaa6 1144 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1145 //
1146 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
1147 allBins[iRow] = new Float_t[maxBin];
adefcaa6 1148 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 1149 allSigBins[iRow] = new Int_t[maxBin];
1150 allNSigBins[iRow]=0;
adefcaa6 1151 }
1152 //
22c352f8 1153 // Loop over sectors
adefcaa6 1154 //
22c352f8 1155 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 1156
22c352f8 1157 Int_t nRows = 0;
1158 Int_t nDDLs = 0, indexDDL = 0;
1159 if (fSector < kNIS) {
1160 nRows = fParam->GetNRowLow();
1161 fSign = (fSector < kNIS/2) ? 1 : -1;
1162 nDDLs = 2;
1163 indexDDL = fSector * 2;
1164 }
1165 else {
1166 nRows = fParam->GetNRowUp();
1167 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1168 nDDLs = 4;
1169 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1170 }
1171
98ee6d31 1172 // load the raw data for corresponding DDLs
1173 rawReader->Reset();
1174 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1175
1176 // select only good sector
1177 input.Next();
1178 if(input.GetSector() != fSector) continue;
1179
1180 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
98ee6d31 1181
22c352f8 1182 for (Int_t iRow = 0; iRow < nRows; iRow++) {
1183 Int_t maxPad;
1184 if (fSector < kNIS)
7c9528d8 1185 maxPad = fParam->GetNPadsLow(iRow);
22c352f8 1186 else
7c9528d8 1187 maxPad = fParam->GetNPadsUp(iRow);
22c352f8 1188
1189 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
22c352f8 1190 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 1191 allNSigBins[iRow] = 0;
22c352f8 1192 }
1193
adefcaa6 1194 Int_t digCounter=0;
22c352f8 1195 // Begin loop over altro data
adefcaa6 1196 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1197 Float_t gain =1;
1198 Int_t lastPad=-1;
98ee6d31 1199
1200 input.Reset();
22c352f8 1201 while (input.Next()) {
22c352f8 1202 if (input.GetSector() != fSector)
7c9528d8 1203 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1204
22c352f8 1205
22c352f8 1206 Int_t iRow = input.GetRow();
9546a7ff 1207 if (iRow < 0){
7c9528d8 1208 continue;
9546a7ff 1209 }
1210
3f0af4ba 1211 if (iRow < 0 || iRow >= nRows){
7c9528d8 1212 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 1213 iRow, 0, nRows -1));
7c9528d8 1214 continue;
3f0af4ba 1215 }
adefcaa6 1216 //pad
1217 Int_t iPad = input.GetPad();
3f0af4ba 1218 if (iPad < 0 || iPad >= nPadsMax) {
7c9528d8 1219 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 1220 iPad, 0, nPadsMax-1));
7c9528d8 1221 continue;
3f0af4ba 1222 }
adefcaa6 1223 if (iPad!=lastPad){
7c9528d8 1224 gain = gainROC->GetValue(iRow,iPad);
1225 lastPad = iPad;
adefcaa6 1226 }
1227 iPad+=3;
1228 //time
1229 Int_t iTimeBin = input.GetTime();
daac2a58 1230 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
7c9528d8 1231 continue;
1232 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 1233 iTimeBin, 0, iTimeBin -1));
daac2a58 1234 }
adefcaa6 1235 iTimeBin+=3;
3f0af4ba 1236
adefcaa6 1237 //signal
22c352f8 1238 Float_t signal = input.GetSignal();
7c9528d8 1239 if (!calcPedestal && signal <= zeroSup) continue;
9546a7ff 1240
940ed1f0 1241 if (!calcPedestal) {
7c9528d8 1242 Int_t bin = iPad*fMaxTime+iTimeBin;
1243 if (gain>0){
1244 allBins[iRow][bin] = signal/gain;
1245 }else{
1246 allBins[iRow][bin] =0;
1247 }
1248 allSigBins[iRow][allNSigBins[iRow]++] = bin;
940ed1f0 1249 }else{
7c9528d8 1250 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
940ed1f0 1251 }
aba7fdfc 1252 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 1253
1254 // Temporary
1255 digCounter++;
22c352f8 1256 } // End of the loop over altro data
adefcaa6 1257 //
1258 //
aba7fdfc 1259 //
1260 //
22c352f8 1261 // Now loop over rows and perform pedestal subtraction
adefcaa6 1262 if (digCounter==0) continue;
7c9528d8 1263 ProcessSectorData(allBins, allSigBins, allNSigBins);
22c352f8 1264 } // End of loop over sectors
98ee6d31 1265
adefcaa6 1266 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1267 delete [] allBins[iRow];
5b33a7f5 1268 delete [] allSigBins[iRow];
adefcaa6 1269 }
1270 delete [] allBins;
5b33a7f5 1271 delete [] allSigBins;
1272 delete [] allNSigBins;
adefcaa6 1273
9546a7ff 1274 if (rawReader->GetEventId() && fOutput ){
98ee6d31 1275 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
9546a7ff 1276 }
ebe95b38 1277
db3576a7 1278 if(rawReader->GetEventId()) {
98ee6d31 1279 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1280 }
1281
1282 if(fBClonesArray) {
1283 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
db3576a7 1284 }
f8aae377 1285}
1286
940ed1f0 1287void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 1288{
940ed1f0 1289
1290 //
1291 // add virtual charge at the edge
1292 //
7041b196 1293 Double_t kMaxDumpSize = 500000;
ebe95b38 1294 if (!fOutput) {
1295 fBDumpSignal =kFALSE;
1296 }else{
1297 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1298 }
1299
f8aae377 1300 fNcluster=0;
f8aae377 1301 fLoop=1;
a1ec4d07 1302 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 1303 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1304 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1305 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1306 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1307 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1308 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 1309 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1310 Int_t i = fSigBins[iSig];
1bfaa9e9 1311 if (i%fMaxTime<=crtime) continue;
5b33a7f5 1312 Float_t *b = &fBins[i];
940ed1f0 1313 //absolute custs
03fe1804 1314 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1315 //
1316 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
aba7fdfc 1317 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1318 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 1319 //
1320 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1321 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1322 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1323 //
1324 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1325 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1326 // sigma cuts
1327 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1328 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1329 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1330
d101caf3 1331 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1332 Int_t dummy=0;
1333 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1334
f8aae377 1335 //}
1336 }
8569a2b0 1337}
65c045f0 1338
eea36fac 1339Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1340 //
679026e6 1341 // Currently hack to filter digital noise (15.06.2008)
eea36fac 1342 // To be parameterized in the AliTPCrecoParam
1343 // More inteligent way to be used in future
679026e6 1344 // Acces to the proper pedestal file needed
eea36fac 1345 //
1346 if (cl->GetMax()<400) return kTRUE;
1347 Double_t ratio = cl->GetQ()/cl->GetMax();
679026e6 1348 if (cl->GetMax()>700){
1349 if ((ratio - int(ratio)>0.8)) return kFALSE;
1350 }
eea36fac 1351 if ((ratio - int(ratio)<0.95)) return kTRUE;
1352 return kFALSE;
eea36fac 1353}
1354
65c045f0 1355
940ed1f0 1356Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1357 //
1358 // process signal on given pad - + streaming of additional information in special mode
1359 //
1360 // id[0] - sector
1361 // id[1] - row
adefcaa6 1362 // id[2] - pad
1363
1364 //
1365 // ESTIMATE pedestal and the noise
1366 //
adefcaa6 1367 const Int_t kPedMax = 100;
1368 Float_t max = 0;
1369 Float_t maxPos = 0;
1370 Int_t median = -1;
1371 Int_t count0 = 0;
1372 Int_t count1 = 0;
940ed1f0 1373 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1374 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
b7a563f9 1375 Int_t firstBin = fRecoParam->GetFirstBin();
65c045f0 1376 //
adefcaa6 1377 UShort_t histo[kPedMax];
f174362f 1378 //memset(histo,0,kPedMax*sizeof(UShort_t));
b842afde 1379 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
adefcaa6 1380 for (Int_t i=0; i<fMaxTime; i++){
1381 if (signal[i]<=0) continue;
940ed1f0 1382 if (signal[i]>max && i>firstBin) {
65c045f0 1383 max = signal[i];
1384 maxPos = i;
adefcaa6 1385 }
1386 if (signal[i]>kPedMax-1) continue;
1387 histo[int(signal[i]+0.5)]++;
1388 count0++;
65c045f0 1389 }
7fef31a6 1390 //
adefcaa6 1391 for (Int_t i=1; i<kPedMax; i++){
1392 if (count1<count0*0.5) median=i;
1393 count1+=histo[i];
1394 }
1395 // truncated mean
65c045f0 1396 //
7041b196 1397 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1398 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1399 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1400 //
1401 for (Int_t idelta=1; idelta<10; idelta++){
1402 if (median-idelta<=0) continue;
1403 if (median+idelta>kPedMax) continue;
1404 if (count06<0.6*count1){
1405 count06+=histo[median-idelta];
1406 mean06 +=histo[median-idelta]*(median-idelta);
1407 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1408 count06+=histo[median+idelta];
1409 mean06 +=histo[median+idelta]*(median+idelta);
1410 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1411 }
1412 if (count09<0.9*count1){
1413 count09+=histo[median-idelta];
1414 mean09 +=histo[median-idelta]*(median-idelta);
1415 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1416 count09+=histo[median+idelta];
1417 mean09 +=histo[median+idelta]*(median+idelta);
1418 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1419 }
adefcaa6 1420 if (count10<0.95*count1){
1421 count10+=histo[median-idelta];
1422 mean +=histo[median-idelta]*(median-idelta);
1423 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1424 count10+=histo[median+idelta];
1425 mean +=histo[median+idelta]*(median+idelta);
1426 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1427 }
1428 }
3f0af4ba 1429 if (count10) {
1430 mean /=count10;
1431 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1432 }
1433 if (count06) {
1434 mean06/=count06;
1435 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1436 }
1437 if (count09) {
1438 mean09/=count09;
1439 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1440 }
940ed1f0 1441 rmsEvent = rms09;
adefcaa6 1442 //
940ed1f0 1443 pedestalEvent = median;
adefcaa6 1444 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1445 //
1446 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1447 //
1448 // Dump mean signal info
1449 //
b2426624 1450 if (AliTPCReconstructor::StreamLevel()>0) {
1451 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1452 "TimeStamp="<<fTimeStamp<<
1453 "EventType="<<fEventType<<
adefcaa6 1454 "Sector="<<uid[0]<<
1455 "Row="<<uid[1]<<
1456 "Pad="<<uid[2]<<
1457 "Max="<<max<<
1458 "MaxPos="<<maxPos<<
1459 //
1460 "Median="<<median<<
1461 "Mean="<<mean<<
1462 "RMS="<<rms<<
1463 "Mean06="<<mean06<<
1464 "RMS06="<<rms06<<
1465 "Mean09="<<mean09<<
1466 "RMS09="<<rms09<<
940ed1f0 1467 "RMSCalib="<<rmsCalib<<
1468 "PedCalib="<<pedestalCalib<<
adefcaa6 1469 "\n";
b2426624 1470 }
adefcaa6 1471 //
1472 // fill pedestal histogram
1473 //
65c045f0 1474 //
adefcaa6 1475 //
1476 //
1477 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1478 Float_t *dsignal = new Float_t[nchannels];
1479 Float_t *dtime = new Float_t[nchannels];
1480 for (Int_t i=0; i<nchannels; i++){
1481 dtime[i] = i;
1482 dsignal[i] = signal[i];
1483 }
03fe1804 1484
11441748 1485 TGraph * graph=0;
7fef31a6 1486 //
a525cc34 1487 // Big signals dumping
1488 //
b2426624 1489 if (AliTPCReconstructor::StreamLevel()>0) {
b7a563f9 1490 if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin())
a525cc34 1491 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1492 "TimeStamp="<<fTimeStamp<<
1493 "EventType="<<fEventType<<
1494 "Sector="<<uid[0]<<
1495 "Row="<<uid[1]<<
1496 "Pad="<<uid[2]<<
1497 "Graph="<<graph<<
1498 "Max="<<max<<
1499 "MaxPos="<<maxPos<<
1500 //
1501 "Median="<<median<<
1502 "Mean="<<mean<<
1503 "RMS="<<rms<<
1504 "Mean06="<<mean06<<
1505 "RMS06="<<rms06<<
1506 "Mean09="<<mean09<<
1507 "RMS09="<<rms09<<
1508 "\n";
1509 delete graph;
b2426624 1510 }
7fef31a6 1511
65c045f0 1512 delete [] dsignal;
1513 delete [] dtime;
940ed1f0 1514 if (rms06>fRecoParam->GetMaxNoise()) {
1515 pedestalEvent+=1024.;
1516 return 1024+median; // sign noisy channel in debug mode
1517 }
65c045f0 1518 return median;
1519}
1520
1521
1522
03fe1804 1523
03fe1804 1524