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