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