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