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