]> 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
Correct sign for calculated b_y (outside measured region).
[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);
7041b196 569 if (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 }
1c53abe2 579
580 fNcluster++;
581}
582
583
584//_____________________________________________________________________________
f8aae377 585void AliTPCclustererMI::Digits2Clusters()
1c53abe2 586{
587 //-----------------------------------------------------------------
588 // This is a simple cluster finder.
589 //-----------------------------------------------------------------
1c53abe2 590
f8aae377 591 if (!fInput) {
592 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 593 return;
594 }
595
f8aae377 596 if (!fOutput) {
597 Error("Digits2Clusters", "output tree not initialised");
598 return;
1c53abe2 599 }
600
13116aec 601 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 602 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
13116aec 603
1c53abe2 604 AliSimDigits digarr, *dummy=&digarr;
605 fRowDig = dummy;
606 fInput->GetBranch("Segment")->SetAddress(&dummy);
607 Stat_t nentries = fInput->GetEntries();
608
f8aae377 609 fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 610
1c53abe2 611 Int_t nclusters = 0;
13116aec 612
1c53abe2 613 for (Int_t n=0; n<nentries; n++) {
614 fInput->GetEvent(n);
508541c7 615 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 616 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
617 continue;
618 }
508541c7 619 Int_t row = fRow;
13116aec 620 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 621 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 622 //
1c53abe2 623 AliTPCClustersRow *clrow= new AliTPCClustersRow();
624 fRowCl = clrow;
625 clrow->SetClass("AliTPCclusterMI");
626 clrow->SetArray(1);
627
628 clrow->SetID(digarr.GetID());
629 fOutput->GetBranch("Segment")->SetAddress(&clrow);
f8aae377 630 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 631
632
f8aae377 633 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 634 fZWidth = fParam->GetZWidth();
635 if (fSector < kNIS) {
f8aae377 636 fMaxPad = fParam->GetNPadsLow(row);
1c53abe2 637 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 638 fPadLength = fParam->GetPadPitchLength(fSector,row);
639 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 640 } else {
f8aae377 641 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 642 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 643 fPadLength = fParam->GetPadPitchLength(fSector,row);
644 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 645 }
646
647
648 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 649 fBins =new Float_t[fMaxBin];
5b33a7f5 650 fSigBins =new Int_t[fMaxBin];
651 fNSigBins = 0;
13116aec 652 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 653
654 if (digarr.First()) //MI change
655 do {
13116aec 656 Float_t dig=digarr.CurrentDigit();
f8aae377 657 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 658 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 659 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 660 Int_t bin = i*fMaxTime+j;
661 fBins[bin]=dig/gain;
662 fSigBins[fNSigBins++]=bin;
1c53abe2 663 } while (digarr.Next());
664 digarr.ExpandTrackBuffer();
665
940ed1f0 666 FindClusters(noiseROC);
8569a2b0 667
668 fOutput->Fill();
88cb7938 669 delete clrow;
670 nclusters+=fNcluster;
5b33a7f5 671 delete[] fBins;
672 delete[] fSigBins;
88cb7938 673 }
f8aae377 674
19dd5b2f 675 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
f8aae377 676}
677
678void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
679{
680//-----------------------------------------------------------------
22c352f8 681// This is a cluster finder for the TPC raw data.
682// The method assumes NO ordering of the altro channels.
683// The pedestal subtraction can be switched on and off
684// using an option of the TPC reconstructor
f8aae377 685//-----------------------------------------------------------------
686
687 if (!fOutput) {
688 Error("Digits2Clusters", "output tree not initialised");
689 return;
690 }
691
f8aae377 692 fRowDig = NULL;
adefcaa6 693 AliTPCROC * roc = AliTPCROC::Instance();
22c352f8 694 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 695 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
696 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
adefcaa6 697 AliTPCRawStream input(rawReader);
ef4ad662 698 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
699 if (fEventHeader){
700 fTimeStamp = fEventHeader->Get("Timestamp");
701 fEventType = fEventHeader->Get("Type");
702 }
703
22c352f8 704
f8aae377 705 Int_t nclusters = 0;
88cb7938 706
f8aae377 707 fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after
708 const Int_t kNIS = fParam->GetNInnerSector();
709 const Int_t kNOS = fParam->GetNOuterSector();
710 const Int_t kNS = kNIS + kNOS;
711 fZWidth = fParam->GetZWidth();
712 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 713 //
714 //alocate memory for sector - maximal case
715 //
22c352f8 716 Float_t** allBins = NULL;
5b33a7f5 717 Int_t** allSigBins = NULL;
718 Int_t* allNSigBins = NULL;
adefcaa6 719 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
720 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
721 allBins = new Float_t*[nRowsMax];
5b33a7f5 722 allSigBins = new Int_t*[nRowsMax];
723 allNSigBins = new Int_t[nRowsMax];
adefcaa6 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];
adefcaa6 728 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 729 allSigBins[iRow] = new Int_t[maxBin];
730 allNSigBins[iRow]=0;
adefcaa6 731 }
732 //
22c352f8 733 // Loop over sectors
adefcaa6 734 //
22c352f8 735 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 736
940ed1f0 737 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
738 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
739 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
22c352f8 740
741 Int_t nRows = 0;
742 Int_t nDDLs = 0, indexDDL = 0;
743 if (fSector < kNIS) {
744 nRows = fParam->GetNRowLow();
745 fSign = (fSector < kNIS/2) ? 1 : -1;
746 nDDLs = 2;
747 indexDDL = fSector * 2;
748 }
749 else {
750 nRows = fParam->GetNRowUp();
751 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
752 nDDLs = 4;
753 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
754 }
755
22c352f8 756 for (Int_t iRow = 0; iRow < nRows; iRow++) {
757 Int_t maxPad;
758 if (fSector < kNIS)
759 maxPad = fParam->GetNPadsLow(iRow);
760 else
761 maxPad = fParam->GetNPadsUp(iRow);
762
763 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
22c352f8 764 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 765 allNSigBins[iRow] = 0;
22c352f8 766 }
767
768 // Loas the raw data for corresponding DDLs
769 rawReader->Reset();
22c352f8 770 input.SetOldRCUFormat(fIsOldRCUFormat);
362c9d61 771 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
adefcaa6 772 Int_t digCounter=0;
22c352f8 773 // Begin loop over altro data
adefcaa6 774 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
775 Float_t gain =1;
776 Int_t lastPad=-1;
22c352f8 777 while (input.Next()) {
adefcaa6 778 digCounter++;
22c352f8 779 if (input.GetSector() != fSector)
780 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
781
22c352f8 782
783 Int_t iRow = input.GetRow();
784 if (iRow < 0 || iRow >= nRows)
785 AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
786 iRow, 0, nRows -1));
adefcaa6 787 //pad
788 Int_t iPad = input.GetPad();
789 if (iPad < 0 || iPad >= nPadsMax)
22c352f8 790 AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 791 iPad, 0, nPadsMax-1));
792 if (iPad!=lastPad){
793 gain = gainROC->GetValue(iRow,iPad);
794 lastPad = iPad;
795 }
796 iPad+=3;
797 //time
798 Int_t iTimeBin = input.GetTime();
799 if ( iTimeBin < 0 || iTimeBin >= fParam->GetMaxTBin())
22c352f8 800 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 801 iTimeBin, 0, iTimeBin -1));
802 iTimeBin+=3;
803 //signal
22c352f8 804 Float_t signal = input.GetSignal();
adefcaa6 805 if (!calcPedestal && signal <= zeroSup) continue;
940ed1f0 806 if (!calcPedestal) {
5b33a7f5 807 Int_t bin = iPad*fMaxTime+iTimeBin;
808 allBins[iRow][bin] = signal/gain;
809 allSigBins[iRow][allNSigBins[iRow]++] = bin;
940ed1f0 810 }else{
811 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
812 }
813 allBins[iRow][iPad*fMaxTime+0]=1.; // pad with signal
22c352f8 814 } // End of the loop over altro data
adefcaa6 815 //
816 //
22c352f8 817 // Now loop over rows and perform pedestal subtraction
adefcaa6 818 if (digCounter==0) continue;
194b0609 819 // if (fPedSubtraction) {
adefcaa6 820 if (calcPedestal) {
22c352f8 821 for (Int_t iRow = 0; iRow < nRows; iRow++) {
822 Int_t maxPad;
823 if (fSector < kNIS)
824 maxPad = fParam->GetNPadsLow(iRow);
825 else
826 maxPad = fParam->GetNPadsUp(iRow);
827
940ed1f0 828 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
829 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 830 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 831 //Float_t pedestal = TMath::Median(fMaxTime, p);
832 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 833 // calib values
834 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
835 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
836 Double_t rmsEvent = rmsCalib;
837 Double_t pedestalEvent = pedestalCalib;
838 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
839 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
840 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
841
842 //
22c352f8 843 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
5b33a7f5 844 Int_t bin = iPad*fMaxTime+iTimeBin;
845 allBins[iRow][bin] -= pedestalEvent;
7fef31a6 846 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
5b33a7f5 847 allBins[iRow][bin] = 0;
7fef31a6 848 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
5b33a7f5 849 allBins[iRow][bin] = 0;
22c352f8 850 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
5b33a7f5 851 allBins[iRow][bin] = 0;
852 if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
853 allBins[iRow][bin] = 0;
854 if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
22c352f8 855 }
856 }
f8aae377 857 }
22c352f8 858 }
22c352f8 859 // Now loop over rows and find clusters
860 for (fRow = 0; fRow < nRows; fRow++) {
861 fRowCl = new AliTPCClustersRow;
862 fRowCl->SetClass("AliTPCclusterMI");
863 fRowCl->SetArray(1);
864 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
865 fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
866
867 fRx = fParam->GetPadRowRadii(fSector, fRow);
868 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 869 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 870 if (fSector < kNIS)
871 fMaxPad = fParam->GetNPadsLow(fRow);
872 else
873 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 874 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 875
876 fBins = allBins[fRow];
5b33a7f5 877 fSigBins = allSigBins[fRow];
878 fNSigBins = allNSigBins[fRow];
22c352f8 879
940ed1f0 880 FindClusters(noiseROC);
22c352f8 881
882 fOutput->Fill();
883 delete fRowCl;
884 nclusters += fNcluster;
885 } // End of loop to find clusters
22c352f8 886 } // End of loop over sectors
adefcaa6 887
888 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
889 delete [] allBins[iRow];
5b33a7f5 890 delete [] allSigBins[iRow];
adefcaa6 891 }
892 delete [] allBins;
5b33a7f5 893 delete [] allSigBins;
894 delete [] allNSigBins;
adefcaa6 895
940ed1f0 896 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
22c352f8 897
f8aae377 898}
899
940ed1f0 900void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 901{
940ed1f0 902
903 //
904 // add virtual charge at the edge
905 //
7041b196 906 Double_t kMaxDumpSize = 500000;
907 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
908 //
f8aae377 909 fNcluster=0;
f8aae377 910 fLoop=1;
a1ec4d07 911 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 912 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
913 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
914 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
915 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
916 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
917 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 918 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
919 Int_t i = fSigBins[iSig];
920 if (i%fMaxTime<=crtime) continue;
921 Float_t *b = &fBins[i];
940ed1f0 922 //absolute custs
03fe1804 923 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
924 //
925 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
926 // if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
927 //if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
928 //
929 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 930 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 931 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 932 //
933 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
934 // sigma cuts
935 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
936 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
937 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
938
03fe1804 939 AliTPCclusterMI c(kFALSE); // default cosntruction without info
f8aae377 940 Int_t dummy=0;
941 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 942
f8aae377 943 //}
944 }
8569a2b0 945}
65c045f0 946
947
940ed1f0 948Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 949 //
950 // process signal on given pad - + streaming of additional information in special mode
951 //
952 // id[0] - sector
953 // id[1] - row
adefcaa6 954 // id[2] - pad
955
956 //
957 // ESTIMATE pedestal and the noise
958 //
adefcaa6 959 const Int_t kPedMax = 100;
7041b196 960 Double_t kMaxDebugSize = 5000000.;
adefcaa6 961 Float_t max = 0;
962 Float_t maxPos = 0;
963 Int_t median = -1;
964 Int_t count0 = 0;
965 Int_t count1 = 0;
940ed1f0 966 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
967 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
968 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 969 //
adefcaa6 970 UShort_t histo[kPedMax];
971 memset(histo,0,kPedMax*sizeof(UShort_t));
972 for (Int_t i=0; i<fMaxTime; i++){
973 if (signal[i]<=0) continue;
940ed1f0 974 if (signal[i]>max && i>firstBin) {
65c045f0 975 max = signal[i];
976 maxPos = i;
adefcaa6 977 }
978 if (signal[i]>kPedMax-1) continue;
979 histo[int(signal[i]+0.5)]++;
980 count0++;
65c045f0 981 }
7fef31a6 982 //
adefcaa6 983 for (Int_t i=1; i<kPedMax; i++){
984 if (count1<count0*0.5) median=i;
985 count1+=histo[i];
986 }
987 // truncated mean
65c045f0 988 //
7041b196 989 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
990 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
991 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 992 //
993 for (Int_t idelta=1; idelta<10; idelta++){
994 if (median-idelta<=0) continue;
995 if (median+idelta>kPedMax) continue;
996 if (count06<0.6*count1){
997 count06+=histo[median-idelta];
998 mean06 +=histo[median-idelta]*(median-idelta);
999 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1000 count06+=histo[median+idelta];
1001 mean06 +=histo[median+idelta]*(median+idelta);
1002 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1003 }
1004 if (count09<0.9*count1){
1005 count09+=histo[median-idelta];
1006 mean09 +=histo[median-idelta]*(median-idelta);
1007 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1008 count09+=histo[median+idelta];
1009 mean09 +=histo[median+idelta]*(median+idelta);
1010 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1011 }
adefcaa6 1012 if (count10<0.95*count1){
1013 count10+=histo[median-idelta];
1014 mean +=histo[median-idelta]*(median-idelta);
1015 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1016 count10+=histo[median+idelta];
1017 mean +=histo[median+idelta]*(median+idelta);
1018 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1019 }
1020 }
1021 mean /=count10;
1022 mean06/=count06;
1023 mean09/=count09;
1024 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1025 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
940ed1f0 1026 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1027 rmsEvent = rms09;
adefcaa6 1028 //
940ed1f0 1029 pedestalEvent = median;
adefcaa6 1030 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1031 //
1032 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1033 //
1034 // Dump mean signal info
1035 //
1036 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1037 "TimeStamp="<<fTimeStamp<<
1038 "EventType="<<fEventType<<
adefcaa6 1039 "Sector="<<uid[0]<<
1040 "Row="<<uid[1]<<
1041 "Pad="<<uid[2]<<
1042 "Max="<<max<<
1043 "MaxPos="<<maxPos<<
1044 //
1045 "Median="<<median<<
1046 "Mean="<<mean<<
1047 "RMS="<<rms<<
1048 "Mean06="<<mean06<<
1049 "RMS06="<<rms06<<
1050 "Mean09="<<mean09<<
1051 "RMS09="<<rms09<<
940ed1f0 1052 "RMSCalib="<<rmsCalib<<
1053 "PedCalib="<<pedestalCalib<<
adefcaa6 1054 "\n";
1055 //
1056 // fill pedestal histogram
1057 //
1058 AliTPCROC * roc = AliTPCROC::Instance();
1059 if (!fAmplitudeHisto){
1060 fAmplitudeHisto = new TObjArray(72);
1061 }
1062 //
1063 if (uid[0]<roc->GetNSectors()
1064 && uid[1]< roc->GetNRows(uid[0]) &&
1065 uid[2] <roc->GetNPads(uid[0], uid[1])){
65c045f0 1066 TObjArray * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
1067 if (!sectorArray){
adefcaa6 1068 Int_t npads =roc->GetNChannels(uid[0]);
65c045f0 1069 sectorArray = new TObjArray(npads);
1070 fAmplitudeHisto->AddAt(sectorArray, uid[0]);
1071 }
adefcaa6 1072 Int_t position = uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
65c045f0 1073 TH1F * histo = (TH1F*)sectorArray->UncheckedAt(position);
1074 if (!histo){
1075 char hname[100];
1076 sprintf(hname,"Amp_%d_%d_%d",uid[0],uid[1],uid[2]);
1077 TFile * backup = gFile;
1078 fDebugStreamer->GetFile()->cd();
194b0609 1079 histo = new TH1F(hname, hname, 100, 5,100);
65c045f0 1080 //histo->SetDirectory(0); // histogram not connected to directory -(File)
1081 sectorArray->AddAt(histo, position);
1082 if (backup) backup->cd();
1083 }
1084 for (Int_t i=0; i<nchannels; i++){
adefcaa6 1085 histo->Fill(signal[i]);
65c045f0 1086 }
1087 }
1088 //
adefcaa6 1089 //
1090 //
1091 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1092 Float_t *dsignal = new Float_t[nchannels];
1093 Float_t *dtime = new Float_t[nchannels];
1094 for (Int_t i=0; i<nchannels; i++){
1095 dtime[i] = i;
1096 dsignal[i] = signal[i];
1097 }
1098 //
03fe1804 1099 // Digital noise
1100 //
d626381a 1101 // if (max-median>30.*TMath::Max(1.,Double_t(rms06)) && (((*fDebugStreamer)<<"SignalDN").GetSize()<kMaxDebugSize)){
1102// //
1103// //
1104// TGraph * graph =new TGraph(nchannels, dtime, dsignal);
1105// //
1106// //
1107// // jumps left - right
1108// Int_t njumps0=0;
1109// Double_t deltaT0[2000];
1110// Double_t deltaA0[2000];
1111// Int_t lastJump0 = fRecoParam->GetFirstBin();
1112// Int_t njumps1=0;
1113// Double_t deltaT1[2000];
1114// Double_t deltaA1[2000];
1115// Int_t lastJump1 = fRecoParam->GetFirstBin();
1116// Int_t njumps2=0;
1117// Double_t deltaT2[2000];
1118// Double_t deltaA2[2000];
1119// Int_t lastJump2 = fRecoParam->GetFirstBin();
1120
1121// for (Int_t itime=fRecoParam->GetFirstBin()+1; itime<fRecoParam->GetLastBin()-1; itime++){
1122// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1123// TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1124// (dsignal[itime-1]-median<5.*rms06) &&
1125// (dsignal[itime+1]-median<5.*rms06)
1126// ){
1127// deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1];
1128// deltaT0[njumps0] = itime-lastJump0;
1129// lastJump0 = itime;
1130// njumps0++;
1131// }
1132// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1133// (dsignal[itime-1]-median<5.*rms06)
1134// ) {
1135// deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1];
1136// deltaT1[njumps1] = itime-lastJump1;
1137// lastJump1 = itime;
1138// njumps1++;
1139// }
1140// if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1141// (dsignal[itime+1]-median<5.*rms06)
1142// ) {
1143// deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1];
1144// deltaT2[njumps2] = itime-lastJump2;
1145// lastJump2 = itime;
1146// njumps2++;
1147// }
1148// }
1149// //
1150// if (njumps0>0 || njumps1>0 || njumps2>0){
1151// TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0);
1152// TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1);
1153// TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2);
1154// (*fDebugStreamer)<<"SignalDN"<< //digital - noise pads - or random sample of pads
1155// "TimeStamp="<<fTimeStamp<<
1156// "EventType="<<fEventType<<
1157// "Sector="<<uid[0]<<
1158// "Row="<<uid[1]<<
1159// "Pad="<<uid[2]<<
1160// "Graph="<<graph<<
1161// "Max="<<max<<
1162// "MaxPos="<<maxPos<<
1163// "Graph.="<<graph<<
1164// "P0GraphDN0.="<<graphDN0<<
1165// "P1GraphDN1.="<<graphDN1<<
1166// "P2GraphDN2.="<<graphDN2<<
1167// //
1168// "Median="<<median<<
1169// "Mean="<<mean<<
1170// "RMS="<<rms<<
1171// "Mean06="<<mean06<<
1172// "RMS06="<<rms06<<
1173// "Mean09="<<mean09<<
1174// "RMS09="<<rms09<<
1175// "\n";
1176// delete graphDN0;
1177// delete graphDN1;
1178// delete graphDN2;
1179// }
1180// delete graph;
1181// }
03fe1804 1182
1183 //
1184 // NOISE STUDY Fourier transform
adefcaa6 1185 //
65c045f0 1186 TGraph * graph;
03fe1804 1187 Bool_t random = (gRandom->Rndm()<0.0003);
7041b196 1188 if (((*fDebugStreamer)<<"SignalN").GetSize()<kMaxDebugSize)
1189 if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
65c045f0 1190 graph =new TGraph(nchannels, dtime, dsignal);
03fe1804 1191 if (rms06>1.*fParam->GetZeroSup() || random){
1192 //Double_t *input, Double_t threshold, Bool_t locMax, Double_t *freq, Double_t *re, Double_t *im, Double_t *mag, Double_t *phi);
1193 Float_t * input = &(dsignal[fRecoParam->GetFirstBin()]);
1194 Float_t freq[2000], re[2000], im[2000], mag[2000], phi[2000];
1195 Int_t npoints = TransformFFT(input, -1,kFALSE, freq, re, im, mag, phi);
1196 TGraph *graphMag0 = new TGraph(npoints, freq, mag);
1197 TGraph *graphPhi0 = new TGraph(npoints, freq, phi);
1198 npoints = TransformFFT(input, 0.5,kTRUE, freq, re, im, mag, phi);
1199 TGraph *graphMag1 = new TGraph(npoints, freq, mag);
1200 TGraph *graphPhi1 = new TGraph(npoints, freq, phi);
1201
65c045f0 1202 (*fDebugStreamer)<<"SignalN"<< //noise pads - or random sample of pads
ef4ad662 1203 "TimeStamp="<<fTimeStamp<<
1204 "EventType="<<fEventType<<
65c045f0 1205 "Sector="<<uid[0]<<
1206 "Row="<<uid[1]<<
1207 "Pad="<<uid[2]<<
03fe1804 1208 "Graph.="<<graph<<
65c045f0 1209 "Max="<<max<<
1210 "MaxPos="<<maxPos<<
1211 //
1212 "Median="<<median<<
1213 "Mean="<<mean<<
1214 "RMS="<<rms<<
1215 "Mean06="<<mean06<<
1216 "RMS06="<<rms06<<
65c045f0 1217 "Mean09="<<mean09<<
1218 "RMS09="<<rms09<<
03fe1804 1219 // FFT part
1220 "Mag0.="<<graphMag0<<
1221 "Mag1.="<<graphMag1<<
1222 "Phi0.="<<graphPhi0<<
1223 "Phi1.="<<graphPhi1<<
65c045f0 1224 "\n";
03fe1804 1225 delete graphMag0;
1226 delete graphMag1;
1227 delete graphPhi0;
1228 delete graphPhi1;
1229 }
1230 //
1231 // Big signals dumping
1232 //
1233
fe5055e5 1234 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
65c045f0 1235 (*fDebugStreamer)<<"SignalB"<< // pads with signal
ef4ad662 1236 "TimeStamp="<<fTimeStamp<<
1237 "EventType="<<fEventType<<
65c045f0 1238 "Sector="<<uid[0]<<
1239 "Row="<<uid[1]<<
1240 "Pad="<<uid[2]<<
864c1803 1241 "Graph="<<graph<<
65c045f0 1242 "Max="<<max<<
1243 "MaxPos="<<maxPos<<
1244 //
1245 "Median="<<median<<
1246 "Mean="<<mean<<
1247 "RMS="<<rms<<
1248 "Mean06="<<mean06<<
1249 "RMS06="<<rms06<<
65c045f0 1250 "Mean09="<<mean09<<
1251 "RMS09="<<rms09<<
1252 "\n";
1253 delete graph;
1254 }
1255
7fef31a6 1256 //
1257 //
1258 // Central Electrode signal analysis
1259 //
7041b196 1260 Float_t ceQmax =0, ceQsum=0, ceTime=0;
1261 Float_t cemean = mean06, cerms=rms06 ;
7fef31a6 1262 Int_t cemaxpos= 0;
7041b196 1263 Float_t ceThreshold=5.*cerms;
1264 Float_t ceSumThreshold=8.*cerms;
97f127bb 1265 const Int_t kCemin=5; // range for the analysis of the ce signal +- channels from the peak
1266 const Int_t kCemax=5;
adefcaa6 1267 for (Int_t i=nchannels-2; i>nchannels/2; i--){
7fef31a6 1268 if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){
1269 cemaxpos=i;
1270 break;
1271 }
1272 }
1273 if (cemaxpos!=0){
adefcaa6 1274 ceQmax = 0;
1275 Int_t cemaxpos2=0;
1276 for (Int_t i=cemaxpos-20; i<cemaxpos+5; i++){
1277 if (i<0 || i>nchannels-1) continue;
1278 Double_t val=dsignal[i]- cemean;
1279 if (val>ceQmax){
1280 cemaxpos2=i;
1281 ceQmax = val;
7fef31a6 1282 }
adefcaa6 1283 }
1284 cemaxpos = cemaxpos2;
1285
1286 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
1287 if (i>0 && i<nchannels&&dsignal[i]- cemean>0){
1288 Double_t val=dsignal[i]- cemean;
1289 ceTime+=val*dtime[i];
1290 ceQsum+=val;
1291 if (val>ceQmax) ceQmax=val;
7fef31a6 1292 }
adefcaa6 1293 }
1294 if (ceQmax&&ceQsum>ceSumThreshold) {
1295 ceTime/=ceQsum;
1296 (*fDebugStreamer)<<"Signalce"<<
ef4ad662 1297 "TimeStamp="<<fTimeStamp<<
1298 "EventType="<<fEventType<<
adefcaa6 1299 "Sector="<<uid[0]<<
1300 "Row="<<uid[1]<<
1301 "Pad="<<uid[2]<<
1302 "Max="<<ceQmax<<
1303 "Qsum="<<ceQsum<<
1304 "Time="<<ceTime<<
1305 "RMS06="<<rms06<<
1306 //
1307 "\n";
1308 }
7fef31a6 1309 }
1310 // end of ce signal analysis
1311 //
1312
1313 //
1314 // Gating grid signal analysis
1315 //
1316 Double_t ggQmax =0, ggQsum=0, ggTime=0;
1317 Double_t ggmean = mean06, ggrms=rms06 ;
1318 Int_t ggmaxpos= 0;
1319 Double_t ggThreshold=5.*ggrms;
1320 Double_t ggSumThreshold=8.*ggrms;
bd85a5e8 1321
adefcaa6 1322 for (Int_t i=1; i<nchannels/4; i++){
bd85a5e8 1323 if ( (dsignal[i]-mean06)>ggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] &&
1324 (dsignal[i]+dsignal[i+1]+dsignal[i-1]-3*mean06)>ggSumThreshold){
7fef31a6 1325 ggmaxpos=i;
bd85a5e8 1326 if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
7fef31a6 1327 break;
1328 }
1329 }
1330 if (ggmaxpos!=0){
bd85a5e8 1331 for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){
1332 if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
7fef31a6 1333 Double_t val=dsignal[i]- ggmean;
1334 ggTime+=val*dtime[i];
1335 ggQsum+=val;
1336 if (val>ggQmax) ggQmax=val;
1337 }
1338 }
1339 if (ggQmax&&ggQsum>ggSumThreshold) {
1340 ggTime/=ggQsum;
1341 (*fDebugStreamer)<<"Signalgg"<<
ef4ad662 1342 "TimeStamp="<<fTimeStamp<<
1343 "EventType="<<fEventType<<
7fef31a6 1344 "Sector="<<uid[0]<<
1345 "Row="<<uid[1]<<
1346 "Pad="<<uid[2]<<
1347 "Max="<<ggQmax<<
1348 "Qsum="<<ggQsum<<
1349 "Time="<<ggTime<<
bd85a5e8 1350 "RMS06="<<rms06<<
7fef31a6 1351 //
1352 "\n";
1353 }
1354 }
1355 // end of gg signal analysis
1356
1357
65c045f0 1358 delete [] dsignal;
1359 delete [] dtime;
940ed1f0 1360 if (rms06>fRecoParam->GetMaxNoise()) {
1361 pedestalEvent+=1024.;
1362 return 1024+median; // sign noisy channel in debug mode
1363 }
65c045f0 1364 return median;
1365}
1366
1367
1368
1369void AliTPCclustererMI::DumpHistos(){
adefcaa6 1370 //
1371 // Dump histogram information
1372 //
65c045f0 1373 if (!fAmplitudeHisto) return;
adefcaa6 1374 AliTPCROC * roc = AliTPCROC::Instance();
65c045f0 1375 for (UInt_t isector=0; isector<AliTPCROC::Instance()->GetNSectors(); isector++){
1376 TObjArray * array = (TObjArray*)fAmplitudeHisto->UncheckedAt(isector);
1377 if (!array) continue;
1378 for (UInt_t ipad = 0; ipad <(UInt_t)array->GetEntriesFast(); ipad++){
1379 TH1F * histo = (TH1F*) array->UncheckedAt(ipad);
1380 if (!histo) continue;
1381 if (histo->GetEntries()<100) continue;
1382 histo->Fit("gaus","q");
1383 Float_t mean = histo->GetMean();
1384 Float_t rms = histo->GetRMS();
1385 Float_t gmean = histo->GetFunction("gaus")->GetParameter(1);
1386 Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2);
03fe1804 1387 Float_t gmeanErr = histo->GetFunction("gaus")->GetParError(1);
1388 Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2);
65c045f0 1389 Float_t max = histo->GetFunction("gaus")->GetParameter(0);
1390
1391 // get pad number
1392 UInt_t row=0, pad =0;
adefcaa6 1393 const UInt_t *indexes =roc->GetRowIndexes(isector);
1394 for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
940ed1f0 1395 if (indexes[irow]<=ipad){
65c045f0 1396 row = irow;
1397 pad = ipad-indexes[irow];
1398 }
1399 }
1400 Int_t rpad = pad - (AliTPCROC::Instance()->GetNPads(isector,row))/2;
1401 //
1402 (*fDebugStreamer)<<"Fit"<<
ef4ad662 1403 "TimeStamp="<<fTimeStamp<<
1404 "EventType="<<fEventType<<
65c045f0 1405 "Sector="<<isector<<
1406 "Row="<<row<<
1407 "Pad="<<pad<<
1408 "RPad="<<rpad<<
1409 "Max="<<max<<
1410 "Mean="<<mean<<
1411 "RMS="<<rms<<
1412 "GMean="<<gmean<<
1413 "GSigma="<<gsigma<<
03fe1804 1414 "GMeanErr="<<gmeanErr<<
1415 "GSigmaErr="<<gsigmaErr<<
65c045f0 1416 "\n";
1417 if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
1418 }
1419 }
1420}
03fe1804 1421
1422
1423
1424Int_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)
1425{
1426 //
1427 // calculate fourrie transform
1428 // return only frequncies with mag over threshold
1429 // if locMax is spectified only freque with local maxima over theshold is returned
1430
1431 if (! fFFTr2c) return kFALSE;
1432 if (!freq) return kFALSE;
1433
1434 Int_t current=0;
1435 Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
1436 Double_t *in = new Double_t[nPoints];
1437 Double_t *rfft = new Double_t[nPoints];
1438 Double_t *ifft = new Double_t[nPoints];
1439 for (Int_t i=0; i<nPoints; i++){in[i]=input[i];}
1440 fFFTr2c->SetPoints(in);
1441 fFFTr2c->Transform();
1442 fFFTr2c->GetPointsComplex(rfft, ifft);
1443 for (Int_t i=3; i<nPoints/2-3; i++){
1444 Float_t lmag = TMath::Sqrt(rfft[i]*rfft[i]+ifft[i]*ifft[i])/nPoints;
1445 if (lmag<threshold) continue;
1446 if (locMax){
1447 if ( TMath::Sqrt(rfft[i-1]*rfft[i-1]+ifft[i-1]*ifft[i-1])/nPoints>lmag) continue;
1448 if ( TMath::Sqrt(rfft[i+1]*rfft[i+1]+ifft[i+1]*ifft[i+1])/nPoints>lmag) continue;
1449 if ( TMath::Sqrt(rfft[i-2]*rfft[i-2]+ifft[i-2]*ifft[i-2])/nPoints>lmag) continue;
1450 if ( TMath::Sqrt(rfft[i+2]*rfft[i+2]+ifft[i+2]*ifft[i+2])/nPoints>lmag) continue;
1451 if ( TMath::Sqrt(rfft[i-3]*rfft[i-3]+ifft[i-3]*ifft[i-3])/nPoints>lmag) continue;
1452 if ( TMath::Sqrt(rfft[i+3]*rfft[i+3]+ifft[i+3]*ifft[i+3])/nPoints>lmag) continue;
1453 }
1454
1455 freq[current] = Float_t(i)/Float_t(nPoints);
1456 //
1457 re[current] = rfft[i];
1458 im[current] = ifft[i];
1459 mag[current]=lmag;
1460 phi[current]=TMath::ATan2(ifft[i],rfft[i]);
1461 current++;
1462 }
1463 delete [] in;
1464 delete [] rfft;
1465 delete [] ifft;
1466 return current;
1467}
1468