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