]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TPC/AliTPCclustererMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Moving Validator to base in order to be able to use it in sim,rec and shuttle
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
... / ...
CommitLineData
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
16/* $Id$ */
17
18//-------------------------------------------------------
19// Implementation of the TPC clusterer
20//
21// Origin: Marian Ivanov
22//-------------------------------------------------------
23
24#include "Riostream.h"
25#include <TF1.h>
26#include <TFile.h>
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>
34
35#include "AliDigits.h"
36#include "AliLoader.h"
37#include "AliLog.h"
38#include "AliMathBase.h"
39#include "AliRawEventHeaderBase.h"
40#include "AliRawReader.h"
41#include "AliRunLoader.h"
42#include "AliSimDigits.h"
43#include "AliTPCCalPad.h"
44#include "AliTPCCalROC.h"
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"
55
56ClassImp(AliTPCclustererMI)
57
58
59
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),
76 fEventHeader(0),
77 fTimeStamp(0),
78 fEventType(0),
79 fInput(0),
80 fOutput(0),
81 fRowCl(0),
82 fRowDig(0),
83 fParam(0),
84 fNcluster(0),
85 fAmplitudeHisto(0),
86 fDebugStreamer(0),
87 fRecoParam(0),
88 fBDumpSignal(kFALSE),
89 fFFTr2c(0)
90{
91 //
92 // COSNTRUCTOR
93 // param - tpc parameters for given file
94 // recoparam - reconstruction parameters
95 //
96 fIsOldRCUFormat = kFALSE;
97 fInput =0;
98 fOutput=0;
99 fParam = par;
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 }
107 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
108 fAmplitudeHisto = 0;
109 Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
110 fFFTr2c = TVirtualFFT::FFT(1, &nPoints, "R2C K");
111}
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//______________________________________________________________
158AliTPCclustererMI::~AliTPCclustererMI(){
159 DumpHistos();
160 if (fAmplitudeHisto) delete fAmplitudeHisto;
161 if (fDebugStreamer) delete fDebugStreamer;
162}
163
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
192 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
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
203 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
204 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
205 Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
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
214void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
215AliTPCclusterMI &c)
216{
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 //
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};
228 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
229 Float_t * resmatrix[5];
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);
241 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
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;
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
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);
322 c.SetPad(meani);
323 c.SetTimeBin(meanj);
324 c.SetSigmaY2(mi2);
325 c.SetSigmaZ2(mj2);
326 AddCluster(c,(Float_t*)vmatrix,k);
327 //remove cluster data from data
328 for (Int_t di=-2;di<=2;di++)
329 for (Int_t dj=-2;dj<=2;dj++){
330 resmatrix[di+2][dj] -= vmatrix[di+2][dj+2];
331 if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
332 }
333 resmatrix[2][0] =0;
334
335 return;
336 }
337 //
338 //unfolding when neccessary
339 //
340
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};
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);
360 c.SetPad(meani);
361 c.SetTimeBin(meanj);
362 c.SetSigmaY2(mi2);
363 c.SetSigmaZ2(mj2);
364 c.SetType(Char_t(overlap)+1);
365 AddCluster(c,(Float_t*)vmatrix,k);
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
376void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
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];
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.);
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];
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.);
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){
486 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
487 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
488 //
489 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
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
523void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
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
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 }
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));
553 if (!fRecoParam->GetBYMirror()){
554 if (fSector%36>17){
555 c.SetY(-(meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
556 }
557 }
558 c.SetZ(fZWidth*(meanj-3));
559 c.SetZ(c.GetZ() - 3.*fParam->GetZSigma() + fParam->GetNTBinsL1()*fParam->GetZWidth()); // PASA delay + L1 delay
560 c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
561 c.SetX(fRx);
562 c.SetDetector(fSector);
563 c.SetRow(fRow);
564
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();
573 AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
574 if (matrix ) {
575 Int_t nbins=0;
576 Float_t *graph =0;
577 if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
578 nbins = fMaxTime;
579 graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
580 }
581 AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
582 cl->SetInfo(info);
583 }
584
585 fNcluster++;
586}
587
588
589//_____________________________________________________________________________
590void AliTPCclustererMI::Digits2Clusters()
591{
592 //-----------------------------------------------------------------
593 // This is a simple cluster finder.
594 //-----------------------------------------------------------------
595
596 if (!fInput) {
597 Error("Digits2Clusters", "input tree not initialised");
598 return;
599 }
600
601 if (!fOutput) {
602 Error("Digits2Clusters", "output tree not initialised");
603 return;
604 }
605
606 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
607 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
608
609 AliSimDigits digarr, *dummy=&digarr;
610 fRowDig = dummy;
611 fInput->GetBranch("Segment")->SetAddress(&dummy);
612 Stat_t nentries = fInput->GetEntries();
613
614 fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after
615
616 Int_t nclusters = 0;
617
618 for (Int_t n=0; n<nentries; n++) {
619 fInput->GetEvent(n);
620 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
621 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
622 continue;
623 }
624 Int_t row = fRow;
625 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
626 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
627 //
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);
635 fRx=fParam->GetPadRowRadii(fSector,row);
636
637
638 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
639 fZWidth = fParam->GetZWidth();
640 if (fSector < kNIS) {
641 fMaxPad = fParam->GetNPadsLow(row);
642 fSign = (fSector < kNIS/2) ? 1 : -1;
643 fPadLength = fParam->GetPadPitchLength(fSector,row);
644 fPadWidth = fParam->GetPadPitchWidth();
645 } else {
646 fMaxPad = fParam->GetNPadsUp(row);
647 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
648 fPadLength = fParam->GetPadPitchLength(fSector,row);
649 fPadWidth = fParam->GetPadPitchWidth();
650 }
651
652
653 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
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);
658
659 if (digarr.First()) //MI change
660 do {
661 Float_t dig=digarr.CurrentDigit();
662 if (dig<=fParam->GetZeroSup()) continue;
663 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
664 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
665 fBins[i*fMaxTime+j]=dig/gain;
666 } while (digarr.Next());
667 digarr.ExpandTrackBuffer();
668
669 FindClusters(noiseROC);
670
671 fOutput->Fill();
672 delete clrow;
673 nclusters+=fNcluster;
674 delete[] fBins;
675 delete[] fResBins;
676 }
677
678 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
679}
680
681void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
682{
683//-----------------------------------------------------------------
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
688//-----------------------------------------------------------------
689
690 if (!fOutput) {
691 Error("Digits2Clusters", "output tree not initialised");
692 return;
693 }
694
695 fRowDig = NULL;
696 AliTPCROC * roc = AliTPCROC::Instance();
697 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
698 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
699 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
700 AliTPCRawStream input(rawReader);
701 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
702 if (fEventHeader){
703 fTimeStamp = fEventHeader->Get("Timestamp");
704 fEventType = fEventHeader->Get("Type");
705 }
706
707
708 Int_t nclusters = 0;
709
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();
716 //
717 //alocate memory for sector - maximal case
718 //
719 Float_t** allBins = NULL;
720 Float_t** allBinsRes = NULL;
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);
731 }
732 //
733 // Loop over sectors
734 //
735 for(fSector = 0; fSector < kNS; fSector++) {
736
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
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
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
764 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
765 }
766
767 // Loas the raw data for corresponding DDLs
768 rawReader->Reset();
769 input.SetOldRCUFormat(fIsOldRCUFormat);
770 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
771 Int_t digCounter=0;
772 // Begin loop over altro data
773 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
774 Float_t gain =1;
775 Int_t lastPad=-1;
776 while (input.Next()) {
777 digCounter++;
778 if (input.GetSector() != fSector)
779 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
780
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));
786 //pad
787 Int_t iPad = input.GetPad();
788 if (iPad < 0 || iPad >= nPadsMax)
789 AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
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())
799 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
800 iTimeBin, 0, iTimeBin -1));
801 iTimeBin+=3;
802 //signal
803 Float_t signal = input.GetSignal();
804 if (!calcPedestal && signal <= zeroSup) continue;
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
811 } // End of the loop over altro data
812 //
813 //
814 // Now loop over rows and perform pedestal subtraction
815 if (digCounter==0) continue;
816 // if (fPedSubtraction) {
817 if (calcPedestal) {
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
825 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
826 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
827 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
828 //Float_t pedestal = TMath::Median(fMaxTime, p);
829 Int_t id[3] = {fSector, iRow, iPad-3};
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 //
840 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
841 allBins[iRow][iPad*fMaxTime+iTimeBin] -= pedestalEvent;
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;
846 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
847 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
848 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < 3.0*rmsEvent) // 3 sigma cut on RMS
849 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
850 }
851 }
852 }
853 }
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);
864 fPadWidth = fParam->GetPadPitchWidth();
865 if (fSector < kNIS)
866 fMaxPad = fParam->GetNPadsLow(fRow);
867 else
868 fMaxPad = fParam->GetNPadsUp(fRow);
869 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
870
871 fBins = allBins[fRow];
872 fResBins = allBinsRes[fRow];
873
874 FindClusters(noiseROC);
875
876 fOutput->Fill();
877 delete fRowCl;
878 nclusters += fNcluster;
879 } // End of loop to find clusters
880 } // End of loop over sectors
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
889 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
890
891}
892
893void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
894{
895
896 //
897 // add virtual charge at the edge
898 //
899 Double_t kMaxDumpSize = 500000;
900 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
901 //
902 if (0) for (Int_t i=0; i<fMaxTime; i++){
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 }
912 fBins[i+2*fMaxTime] = amp0;
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 }
922 fBins[(fMaxPad+3)*fMaxTime+i] = amp0;
923 }
924 memcpy(fResBins,fBins, fMaxBin*sizeof(Float_t));
925 //
926 //
927 //
928 fNcluster=0;
929 fLoop=1;
930 Float_t *b=&fBins[-1]+2*fMaxTime;
931 Int_t crtime = Int_t((fParam->GetZLength()-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
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();
938 for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
939 b++;
940 if (i%fMaxTime<crtime) {
941 Int_t delta = -(i%fMaxTime)+crtime;
942 b+=delta;
943 i+=delta;
944 continue;
945 }
946 //absolute custs
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)
954 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
955 if (!IsMaximum(*b,fMaxTime,b)) continue;
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
963 AliTPCclusterMI c(kFALSE); // default cosntruction without info
964 Int_t dummy=0;
965 MakeCluster(i, fMaxTime, fBins, dummy,c);
966
967 //}
968 }
969}
970
971
972Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
973 //
974 // process signal on given pad - + streaming of additional information in special mode
975 //
976 // id[0] - sector
977 // id[1] - row
978 // id[2] - pad
979
980 //
981 // ESTIMATE pedestal and the noise
982 //
983 const Int_t kPedMax = 100;
984 Double_t kMaxDebugSize = 5000000.;
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;
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();
993 //
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;
998 if (signal[i]>max && i>firstBin) {
999 max = signal[i];
1000 maxPos = i;
1001 }
1002 if (signal[i]>kPedMax-1) continue;
1003 histo[int(signal[i]+0.5)]++;
1004 count0++;
1005 }
1006 //
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
1012 //
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 ;
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);
1035 }
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));
1050 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1051 rmsEvent = rms09;
1052 //
1053 pedestalEvent = median;
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"<<
1061 "TimeStamp="<<fTimeStamp<<
1062 "EventType="<<fEventType<<
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<<
1076 "RMSCalib="<<rmsCalib<<
1077 "PedCalib="<<pedestalCalib<<
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])){
1090 TObjArray * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
1091 if (!sectorArray){
1092 Int_t npads =roc->GetNChannels(uid[0]);
1093 sectorArray = new TObjArray(npads);
1094 fAmplitudeHisto->AddAt(sectorArray, uid[0]);
1095 }
1096 Int_t position = uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
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();
1103 histo = new TH1F(hname, hname, 100, 5,100);
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++){
1109 histo->Fill(signal[i]);
1110 }
1111 }
1112 //
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 //
1123 // Digital noise
1124 //
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 }
1206
1207 //
1208 // NOISE STUDY Fourier transform
1209 //
1210 TGraph * graph;
1211 Bool_t random = (gRandom->Rndm()<0.0003);
1212 if (((*fDebugStreamer)<<"SignalN").GetSize()<kMaxDebugSize)
1213 if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
1214 graph =new TGraph(nchannels, dtime, dsignal);
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
1226 (*fDebugStreamer)<<"SignalN"<< //noise pads - or random sample of pads
1227 "TimeStamp="<<fTimeStamp<<
1228 "EventType="<<fEventType<<
1229 "Sector="<<uid[0]<<
1230 "Row="<<uid[1]<<
1231 "Pad="<<uid[2]<<
1232 "Graph.="<<graph<<
1233 "Max="<<max<<
1234 "MaxPos="<<maxPos<<
1235 //
1236 "Median="<<median<<
1237 "Mean="<<mean<<
1238 "RMS="<<rms<<
1239 "Mean06="<<mean06<<
1240 "RMS06="<<rms06<<
1241 "Mean09="<<mean09<<
1242 "RMS09="<<rms09<<
1243 // FFT part
1244 "Mag0.="<<graphMag0<<
1245 "Mag1.="<<graphMag1<<
1246 "Phi0.="<<graphPhi0<<
1247 "Phi1.="<<graphPhi1<<
1248 "\n";
1249 delete graphMag0;
1250 delete graphMag1;
1251 delete graphPhi0;
1252 delete graphPhi1;
1253 }
1254 //
1255 // Big signals dumping
1256 //
1257
1258 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1259 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1260 "TimeStamp="<<fTimeStamp<<
1261 "EventType="<<fEventType<<
1262 "Sector="<<uid[0]<<
1263 "Row="<<uid[1]<<
1264 "Pad="<<uid[2]<<
1265 "Graph="<<graph<<
1266 "Max="<<max<<
1267 "MaxPos="<<maxPos<<
1268 //
1269 "Median="<<median<<
1270 "Mean="<<mean<<
1271 "RMS="<<rms<<
1272 "Mean06="<<mean06<<
1273 "RMS06="<<rms06<<
1274 "Mean09="<<mean09<<
1275 "RMS09="<<rms09<<
1276 "\n";
1277 delete graph;
1278 }
1279
1280 //
1281 //
1282 // Central Electrode signal analysis
1283 //
1284 Float_t ceQmax =0, ceQsum=0, ceTime=0;
1285 Float_t cemean = mean06, cerms=rms06 ;
1286 Int_t cemaxpos= 0;
1287 Float_t ceThreshold=5.*cerms;
1288 Float_t ceSumThreshold=8.*cerms;
1289 const Int_t kCemin=5; // range for the analysis of the ce signal +- channels from the peak
1290 const Int_t kCemax=5;
1291 for (Int_t i=nchannels-2; i>nchannels/2; i--){
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){
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;
1306 }
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;
1316 }
1317 }
1318 if (ceQmax&&ceQsum>ceSumThreshold) {
1319 ceTime/=ceQsum;
1320 (*fDebugStreamer)<<"Signalce"<<
1321 "TimeStamp="<<fTimeStamp<<
1322 "EventType="<<fEventType<<
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 }
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;
1345
1346 for (Int_t i=1; i<nchannels/4; i++){
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){
1349 ggmaxpos=i;
1350 if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
1351 break;
1352 }
1353 }
1354 if (ggmaxpos!=0){
1355 for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){
1356 if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
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"<<
1366 "TimeStamp="<<fTimeStamp<<
1367 "EventType="<<fEventType<<
1368 "Sector="<<uid[0]<<
1369 "Row="<<uid[1]<<
1370 "Pad="<<uid[2]<<
1371 "Max="<<ggQmax<<
1372 "Qsum="<<ggQsum<<
1373 "Time="<<ggTime<<
1374 "RMS06="<<rms06<<
1375 //
1376 "\n";
1377 }
1378 }
1379 // end of gg signal analysis
1380
1381
1382 delete [] dsignal;
1383 delete [] dtime;
1384 if (rms06>fRecoParam->GetMaxNoise()) {
1385 pedestalEvent+=1024.;
1386 return 1024+median; // sign noisy channel in debug mode
1387 }
1388 return median;
1389}
1390
1391
1392
1393void AliTPCclustererMI::DumpHistos(){
1394 //
1395 // Dump histogram information
1396 //
1397 if (!fAmplitudeHisto) return;
1398 AliTPCROC * roc = AliTPCROC::Instance();
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);
1411 Float_t gmeanErr = histo->GetFunction("gaus")->GetParError(1);
1412 Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2);
1413 Float_t max = histo->GetFunction("gaus")->GetParameter(0);
1414
1415 // get pad number
1416 UInt_t row=0, pad =0;
1417 const UInt_t *indexes =roc->GetRowIndexes(isector);
1418 for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
1419 if (indexes[irow]<=ipad){
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"<<
1427 "TimeStamp="<<fTimeStamp<<
1428 "EventType="<<fEventType<<
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<<
1438 "GMeanErr="<<gmeanErr<<
1439 "GSigmaErr="<<gsigmaErr<<
1440 "\n";
1441 if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
1442 }
1443 }
1444}
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