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