]> 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
New mapping in agreement with the new instructions from Paolo and Giacinto+tools...
[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 }
aba7fdfc 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 //
aba7fdfc 834 //
835 //
22c352f8 836 // Now loop over rows and perform pedestal subtraction
adefcaa6 837 if (digCounter==0) continue;
aba7fdfc 838 // if (calcPedestal) {
839 if (kTRUE) {
22c352f8 840 for (Int_t iRow = 0; iRow < nRows; iRow++) {
841 Int_t maxPad;
842 if (fSector < kNIS)
843 maxPad = fParam->GetNPadsLow(iRow);
844 else
845 maxPad = fParam->GetNPadsUp(iRow);
846
940ed1f0 847 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
aba7fdfc 848 //
849 // Temporary fix for data production - !!!! MARIAN
850 // The noise calibration should take mean and RMS - currently the Gaussian fit used
851 // In case of double peak - the pad should be rejected
852 //
853 // Line mean - if more than given digits over threshold - make a noise calculation
854 // and pedestal substration
855 if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
856 //
940ed1f0 857 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 858 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 859 //Float_t pedestal = TMath::Median(fMaxTime, p);
860 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 861 // calib values
862 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
863 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
864 Double_t rmsEvent = rmsCalib;
865 Double_t pedestalEvent = pedestalCalib;
866 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
867 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
868 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
869
870 //
22c352f8 871 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
5b33a7f5 872 Int_t bin = iPad*fMaxTime+iTimeBin;
873 allBins[iRow][bin] -= pedestalEvent;
7fef31a6 874 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
5b33a7f5 875 allBins[iRow][bin] = 0;
7fef31a6 876 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
5b33a7f5 877 allBins[iRow][bin] = 0;
22c352f8 878 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
5b33a7f5 879 allBins[iRow][bin] = 0;
880 if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
881 allBins[iRow][bin] = 0;
882 if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
22c352f8 883 }
884 }
f8aae377 885 }
22c352f8 886 }
22c352f8 887 // Now loop over rows and find clusters
888 for (fRow = 0; fRow < nRows; fRow++) {
889 fRowCl = new AliTPCClustersRow;
890 fRowCl->SetClass("AliTPCclusterMI");
891 fRowCl->SetArray(1);
892 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
893 fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
894
895 fRx = fParam->GetPadRowRadii(fSector, fRow);
896 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 897 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 898 if (fSector < kNIS)
899 fMaxPad = fParam->GetNPadsLow(fRow);
900 else
901 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 902 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 903
904 fBins = allBins[fRow];
5b33a7f5 905 fSigBins = allSigBins[fRow];
906 fNSigBins = allNSigBins[fRow];
22c352f8 907
940ed1f0 908 FindClusters(noiseROC);
22c352f8 909
910 fOutput->Fill();
911 delete fRowCl;
912 nclusters += fNcluster;
913 } // End of loop to find clusters
22c352f8 914 } // End of loop over sectors
adefcaa6 915
916 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
917 delete [] allBins[iRow];
5b33a7f5 918 delete [] allSigBins[iRow];
adefcaa6 919 }
920 delete [] allBins;
5b33a7f5 921 delete [] allSigBins;
922 delete [] allNSigBins;
adefcaa6 923
f2b386d5 924 if (rawReader->GetEventId() && fOutput ){
925 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
926 }
22c352f8 927
f8aae377 928}
929
940ed1f0 930void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 931{
940ed1f0 932
933 //
934 // add virtual charge at the edge
935 //
7041b196 936 Double_t kMaxDumpSize = 500000;
937 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
938 //
f8aae377 939 fNcluster=0;
f8aae377 940 fLoop=1;
a1ec4d07 941 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 942 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
943 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
944 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
945 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
946 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
947 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 948 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
949 Int_t i = fSigBins[iSig];
950 if (i%fMaxTime<=crtime) continue;
951 Float_t *b = &fBins[i];
940ed1f0 952 //absolute custs
03fe1804 953 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
954 //
955 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
aba7fdfc 956 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
957 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 958 //
959 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 960 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 961 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 962 //
963 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
964 // sigma cuts
965 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
966 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
967 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
968
03fe1804 969 AliTPCclusterMI c(kFALSE); // default cosntruction without info
f8aae377 970 Int_t dummy=0;
971 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 972
f8aae377 973 //}
974 }
8569a2b0 975}
65c045f0 976
977
940ed1f0 978Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 979 //
980 // process signal on given pad - + streaming of additional information in special mode
981 //
982 // id[0] - sector
983 // id[1] - row
adefcaa6 984 // id[2] - pad
985
986 //
987 // ESTIMATE pedestal and the noise
988 //
adefcaa6 989 const Int_t kPedMax = 100;
7041b196 990 Double_t kMaxDebugSize = 5000000.;
adefcaa6 991 Float_t max = 0;
992 Float_t maxPos = 0;
993 Int_t median = -1;
994 Int_t count0 = 0;
995 Int_t count1 = 0;
940ed1f0 996 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
997 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
998 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 999 //
adefcaa6 1000 UShort_t histo[kPedMax];
1001 memset(histo,0,kPedMax*sizeof(UShort_t));
1002 for (Int_t i=0; i<fMaxTime; i++){
1003 if (signal[i]<=0) continue;
940ed1f0 1004 if (signal[i]>max && i>firstBin) {
65c045f0 1005 max = signal[i];
1006 maxPos = i;
adefcaa6 1007 }
1008 if (signal[i]>kPedMax-1) continue;
1009 histo[int(signal[i]+0.5)]++;
1010 count0++;
65c045f0 1011 }
7fef31a6 1012 //
adefcaa6 1013 for (Int_t i=1; i<kPedMax; i++){
1014 if (count1<count0*0.5) median=i;
1015 count1+=histo[i];
1016 }
1017 // truncated mean
65c045f0 1018 //
7041b196 1019 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1020 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1021 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1022 //
1023 for (Int_t idelta=1; idelta<10; idelta++){
1024 if (median-idelta<=0) continue;
1025 if (median+idelta>kPedMax) continue;
1026 if (count06<0.6*count1){
1027 count06+=histo[median-idelta];
1028 mean06 +=histo[median-idelta]*(median-idelta);
1029 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1030 count06+=histo[median+idelta];
1031 mean06 +=histo[median+idelta]*(median+idelta);
1032 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1033 }
1034 if (count09<0.9*count1){
1035 count09+=histo[median-idelta];
1036 mean09 +=histo[median-idelta]*(median-idelta);
1037 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1038 count09+=histo[median+idelta];
1039 mean09 +=histo[median+idelta]*(median+idelta);
1040 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1041 }
adefcaa6 1042 if (count10<0.95*count1){
1043 count10+=histo[median-idelta];
1044 mean +=histo[median-idelta]*(median-idelta);
1045 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1046 count10+=histo[median+idelta];
1047 mean +=histo[median+idelta]*(median+idelta);
1048 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1049 }
1050 }
3f0af4ba 1051 if (count10) {
1052 mean /=count10;
1053 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1054 }
1055 if (count06) {
1056 mean06/=count06;
1057 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1058 }
1059 if (count09) {
1060 mean09/=count09;
1061 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1062 }
940ed1f0 1063 rmsEvent = rms09;
adefcaa6 1064 //
940ed1f0 1065 pedestalEvent = median;
adefcaa6 1066 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1067 //
1068 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1069 //
1070 // Dump mean signal info
1071 //
1072 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1073 "TimeStamp="<<fTimeStamp<<
1074 "EventType="<<fEventType<<
adefcaa6 1075 "Sector="<<uid[0]<<
1076 "Row="<<uid[1]<<
1077 "Pad="<<uid[2]<<
1078 "Max="<<max<<
1079 "MaxPos="<<maxPos<<
1080 //
1081 "Median="<<median<<
1082 "Mean="<<mean<<
1083 "RMS="<<rms<<
1084 "Mean06="<<mean06<<
1085 "RMS06="<<rms06<<
1086 "Mean09="<<mean09<<
1087 "RMS09="<<rms09<<
940ed1f0 1088 "RMSCalib="<<rmsCalib<<
1089 "PedCalib="<<pedestalCalib<<
adefcaa6 1090 "\n";
1091 //
1092 // fill pedestal histogram
1093 //
1094 AliTPCROC * roc = AliTPCROC::Instance();
1095 if (!fAmplitudeHisto){
1096 fAmplitudeHisto = new TObjArray(72);
1097 }
1098 //
1099 if (uid[0]<roc->GetNSectors()
1100 && uid[1]< roc->GetNRows(uid[0]) &&
1101 uid[2] <roc->GetNPads(uid[0], uid[1])){
65c045f0 1102 TObjArray * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
1103 if (!sectorArray){
adefcaa6 1104 Int_t npads =roc->GetNChannels(uid[0]);
65c045f0 1105 sectorArray = new TObjArray(npads);
1106 fAmplitudeHisto->AddAt(sectorArray, uid[0]);
1107 }
adefcaa6 1108 Int_t position = uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
3ad5a1ce 1109 // TH1F * histo = (TH1F*)sectorArray->UncheckedAt(position);
1110// if (!histo){
1111// char hname[100];
1112// sprintf(hname,"Amp_%d_%d_%d",uid[0],uid[1],uid[2]);
1113// TFile * backup = gFile;
1114// fDebugStreamer->GetFile()->cd();
1115// histo = new TH1F(hname, hname, 100, 5,100);
1116// //histo->SetDirectory(0); // histogram not connected to directory -(File)
1117// sectorArray->AddAt(histo, position);
1118// if (backup) backup->cd();
1119// }
1120// for (Int_t i=0; i<nchannels; i++){
1121// histo->Fill(signal[i]);
1122// }
65c045f0 1123 }
1124 //
adefcaa6 1125 //
1126 //
1127 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1128 Float_t *dsignal = new Float_t[nchannels];
1129 Float_t *dtime = new Float_t[nchannels];
1130 for (Int_t i=0; i<nchannels; i++){
1131 dtime[i] = i;
1132 dsignal[i] = signal[i];
1133 }
1134 //
03fe1804 1135 // Digital noise
1136 //
d626381a 1137 // if (max-median>30.*TMath::Max(1.,Double_t(rms06)) && (((*fDebugStreamer)<<"SignalDN").GetSize()<kMaxDebugSize)){
1138// //
1139// //
1140// TGraph * graph =new TGraph(nchannels, dtime, dsignal);
1141// //
1142// //
1143// // jumps left - right
1144// Int_t njumps0=0;
1145// Double_t deltaT0[2000];
1146// Double_t deltaA0[2000];
1147// Int_t lastJump0 = fRecoParam->GetFirstBin();
1148// Int_t njumps1=0;
1149// Double_t deltaT1[2000];
1150// Double_t deltaA1[2000];
1151// Int_t lastJump1 = fRecoParam->GetFirstBin();
1152// Int_t njumps2=0;
1153// Double_t deltaT2[2000];
1154// Double_t deltaA2[2000];
1155// Int_t lastJump2 = fRecoParam->GetFirstBin();
1156
1157// for (Int_t itime=fRecoParam->GetFirstBin()+1; itime<fRecoParam->GetLastBin()-1; itime++){
1158// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1159// TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1160// (dsignal[itime-1]-median<5.*rms06) &&
1161// (dsignal[itime+1]-median<5.*rms06)
1162// ){
1163// deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1];
1164// deltaT0[njumps0] = itime-lastJump0;
1165// lastJump0 = itime;
1166// njumps0++;
1167// }
1168// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1169// (dsignal[itime-1]-median<5.*rms06)
1170// ) {
1171// deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1];
1172// deltaT1[njumps1] = itime-lastJump1;
1173// lastJump1 = itime;
1174// njumps1++;
1175// }
1176// if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1177// (dsignal[itime+1]-median<5.*rms06)
1178// ) {
1179// deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1];
1180// deltaT2[njumps2] = itime-lastJump2;
1181// lastJump2 = itime;
1182// njumps2++;
1183// }
1184// }
1185// //
1186// if (njumps0>0 || njumps1>0 || njumps2>0){
1187// TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0);
1188// TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1);
1189// TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2);
1190// (*fDebugStreamer)<<"SignalDN"<< //digital - noise pads - or random sample of pads
1191// "TimeStamp="<<fTimeStamp<<
1192// "EventType="<<fEventType<<
1193// "Sector="<<uid[0]<<
1194// "Row="<<uid[1]<<
1195// "Pad="<<uid[2]<<
1196// "Graph="<<graph<<
1197// "Max="<<max<<
1198// "MaxPos="<<maxPos<<
1199// "Graph.="<<graph<<
1200// "P0GraphDN0.="<<graphDN0<<
1201// "P1GraphDN1.="<<graphDN1<<
1202// "P2GraphDN2.="<<graphDN2<<
1203// //
1204// "Median="<<median<<
1205// "Mean="<<mean<<
1206// "RMS="<<rms<<
1207// "Mean06="<<mean06<<
1208// "RMS06="<<rms06<<
1209// "Mean09="<<mean09<<
1210// "RMS09="<<rms09<<
1211// "\n";
1212// delete graphDN0;
1213// delete graphDN1;
1214// delete graphDN2;
1215// }
1216// delete graph;
1217// }
03fe1804 1218
1219 //
1220 // NOISE STUDY Fourier transform
adefcaa6 1221 //
65c045f0 1222 TGraph * graph;
03fe1804 1223 Bool_t random = (gRandom->Rndm()<0.0003);
7041b196 1224 if (((*fDebugStreamer)<<"SignalN").GetSize()<kMaxDebugSize)
1225 if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
65c045f0 1226 graph =new TGraph(nchannels, dtime, dsignal);
03fe1804 1227 if (rms06>1.*fParam->GetZeroSup() || random){
1228 //Double_t *input, Double_t threshold, Bool_t locMax, Double_t *freq, Double_t *re, Double_t *im, Double_t *mag, Double_t *phi);
1229 Float_t * input = &(dsignal[fRecoParam->GetFirstBin()]);
1230 Float_t freq[2000], re[2000], im[2000], mag[2000], phi[2000];
1231 Int_t npoints = TransformFFT(input, -1,kFALSE, freq, re, im, mag, phi);
1232 TGraph *graphMag0 = new TGraph(npoints, freq, mag);
1233 TGraph *graphPhi0 = new TGraph(npoints, freq, phi);
1234 npoints = TransformFFT(input, 0.5,kTRUE, freq, re, im, mag, phi);
1235 TGraph *graphMag1 = new TGraph(npoints, freq, mag);
1236 TGraph *graphPhi1 = new TGraph(npoints, freq, phi);
1237
65c045f0 1238 (*fDebugStreamer)<<"SignalN"<< //noise pads - or random sample of pads
ef4ad662 1239 "TimeStamp="<<fTimeStamp<<
1240 "EventType="<<fEventType<<
65c045f0 1241 "Sector="<<uid[0]<<
1242 "Row="<<uid[1]<<
1243 "Pad="<<uid[2]<<
03fe1804 1244 "Graph.="<<graph<<
65c045f0 1245 "Max="<<max<<
1246 "MaxPos="<<maxPos<<
1247 //
1248 "Median="<<median<<
1249 "Mean="<<mean<<
1250 "RMS="<<rms<<
1251 "Mean06="<<mean06<<
1252 "RMS06="<<rms06<<
65c045f0 1253 "Mean09="<<mean09<<
1254 "RMS09="<<rms09<<
03fe1804 1255 // FFT part
1256 "Mag0.="<<graphMag0<<
1257 "Mag1.="<<graphMag1<<
1258 "Phi0.="<<graphPhi0<<
1259 "Phi1.="<<graphPhi1<<
65c045f0 1260 "\n";
03fe1804 1261 delete graphMag0;
1262 delete graphMag1;
1263 delete graphPhi0;
1264 delete graphPhi1;
1265 }
1266 //
1267 // Big signals dumping
1268 //
1269
fe5055e5 1270 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
65c045f0 1271 (*fDebugStreamer)<<"SignalB"<< // pads with signal
ef4ad662 1272 "TimeStamp="<<fTimeStamp<<
1273 "EventType="<<fEventType<<
65c045f0 1274 "Sector="<<uid[0]<<
1275 "Row="<<uid[1]<<
1276 "Pad="<<uid[2]<<
864c1803 1277 "Graph="<<graph<<
65c045f0 1278 "Max="<<max<<
1279 "MaxPos="<<maxPos<<
1280 //
1281 "Median="<<median<<
1282 "Mean="<<mean<<
1283 "RMS="<<rms<<
1284 "Mean06="<<mean06<<
1285 "RMS06="<<rms06<<
65c045f0 1286 "Mean09="<<mean09<<
1287 "RMS09="<<rms09<<
1288 "\n";
1289 delete graph;
1290 }
1291
7fef31a6 1292 //
1293 //
1294 // Central Electrode signal analysis
1295 //
7041b196 1296 Float_t ceQmax =0, ceQsum=0, ceTime=0;
1297 Float_t cemean = mean06, cerms=rms06 ;
7fef31a6 1298 Int_t cemaxpos= 0;
7041b196 1299 Float_t ceThreshold=5.*cerms;
1300 Float_t ceSumThreshold=8.*cerms;
97f127bb 1301 const Int_t kCemin=5; // range for the analysis of the ce signal +- channels from the peak
1302 const Int_t kCemax=5;
adefcaa6 1303 for (Int_t i=nchannels-2; i>nchannels/2; i--){
7fef31a6 1304 if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){
1305 cemaxpos=i;
1306 break;
1307 }
1308 }
1309 if (cemaxpos!=0){
adefcaa6 1310 ceQmax = 0;
1311 Int_t cemaxpos2=0;
1312 for (Int_t i=cemaxpos-20; i<cemaxpos+5; i++){
1313 if (i<0 || i>nchannels-1) continue;
1314 Double_t val=dsignal[i]- cemean;
1315 if (val>ceQmax){
1316 cemaxpos2=i;
1317 ceQmax = val;
7fef31a6 1318 }
adefcaa6 1319 }
1320 cemaxpos = cemaxpos2;
1321
1322 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
1323 if (i>0 && i<nchannels&&dsignal[i]- cemean>0){
1324 Double_t val=dsignal[i]- cemean;
1325 ceTime+=val*dtime[i];
1326 ceQsum+=val;
1327 if (val>ceQmax) ceQmax=val;
7fef31a6 1328 }
adefcaa6 1329 }
1330 if (ceQmax&&ceQsum>ceSumThreshold) {
1331 ceTime/=ceQsum;
1332 (*fDebugStreamer)<<"Signalce"<<
ef4ad662 1333 "TimeStamp="<<fTimeStamp<<
1334 "EventType="<<fEventType<<
adefcaa6 1335 "Sector="<<uid[0]<<
1336 "Row="<<uid[1]<<
1337 "Pad="<<uid[2]<<
1338 "Max="<<ceQmax<<
1339 "Qsum="<<ceQsum<<
1340 "Time="<<ceTime<<
1341 "RMS06="<<rms06<<
1342 //
1343 "\n";
1344 }
7fef31a6 1345 }
1346 // end of ce signal analysis
1347 //
1348
1349 //
1350 // Gating grid signal analysis
1351 //
1352 Double_t ggQmax =0, ggQsum=0, ggTime=0;
1353 Double_t ggmean = mean06, ggrms=rms06 ;
1354 Int_t ggmaxpos= 0;
1355 Double_t ggThreshold=5.*ggrms;
1356 Double_t ggSumThreshold=8.*ggrms;
bd85a5e8 1357
adefcaa6 1358 for (Int_t i=1; i<nchannels/4; i++){
bd85a5e8 1359 if ( (dsignal[i]-mean06)>ggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] &&
1360 (dsignal[i]+dsignal[i+1]+dsignal[i-1]-3*mean06)>ggSumThreshold){
7fef31a6 1361 ggmaxpos=i;
bd85a5e8 1362 if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
7fef31a6 1363 break;
1364 }
1365 }
1366 if (ggmaxpos!=0){
bd85a5e8 1367 for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){
1368 if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
7fef31a6 1369 Double_t val=dsignal[i]- ggmean;
1370 ggTime+=val*dtime[i];
1371 ggQsum+=val;
1372 if (val>ggQmax) ggQmax=val;
1373 }
1374 }
1375 if (ggQmax&&ggQsum>ggSumThreshold) {
1376 ggTime/=ggQsum;
1377 (*fDebugStreamer)<<"Signalgg"<<
ef4ad662 1378 "TimeStamp="<<fTimeStamp<<
1379 "EventType="<<fEventType<<
7fef31a6 1380 "Sector="<<uid[0]<<
1381 "Row="<<uid[1]<<
1382 "Pad="<<uid[2]<<
1383 "Max="<<ggQmax<<
1384 "Qsum="<<ggQsum<<
1385 "Time="<<ggTime<<
bd85a5e8 1386 "RMS06="<<rms06<<
7fef31a6 1387 //
1388 "\n";
1389 }
1390 }
1391 // end of gg signal analysis
1392
1393
65c045f0 1394 delete [] dsignal;
1395 delete [] dtime;
940ed1f0 1396 if (rms06>fRecoParam->GetMaxNoise()) {
1397 pedestalEvent+=1024.;
1398 return 1024+median; // sign noisy channel in debug mode
1399 }
65c045f0 1400 return median;
1401}
1402
1403
1404
1405void AliTPCclustererMI::DumpHistos(){
adefcaa6 1406 //
1407 // Dump histogram information
1408 //
65c045f0 1409 if (!fAmplitudeHisto) return;
adefcaa6 1410 AliTPCROC * roc = AliTPCROC::Instance();
65c045f0 1411 for (UInt_t isector=0; isector<AliTPCROC::Instance()->GetNSectors(); isector++){
1412 TObjArray * array = (TObjArray*)fAmplitudeHisto->UncheckedAt(isector);
1413 if (!array) continue;
1414 for (UInt_t ipad = 0; ipad <(UInt_t)array->GetEntriesFast(); ipad++){
1415 TH1F * histo = (TH1F*) array->UncheckedAt(ipad);
1416 if (!histo) continue;
1417 if (histo->GetEntries()<100) continue;
1418 histo->Fit("gaus","q");
1419 Float_t mean = histo->GetMean();
1420 Float_t rms = histo->GetRMS();
1421 Float_t gmean = histo->GetFunction("gaus")->GetParameter(1);
1422 Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2);
03fe1804 1423 Float_t gmeanErr = histo->GetFunction("gaus")->GetParError(1);
1424 Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2);
65c045f0 1425 Float_t max = histo->GetFunction("gaus")->GetParameter(0);
1426
1427 // get pad number
1428 UInt_t row=0, pad =0;
adefcaa6 1429 const UInt_t *indexes =roc->GetRowIndexes(isector);
1430 for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
940ed1f0 1431 if (indexes[irow]<=ipad){
65c045f0 1432 row = irow;
1433 pad = ipad-indexes[irow];
1434 }
1435 }
1436 Int_t rpad = pad - (AliTPCROC::Instance()->GetNPads(isector,row))/2;
1437 //
1438 (*fDebugStreamer)<<"Fit"<<
ef4ad662 1439 "TimeStamp="<<fTimeStamp<<
1440 "EventType="<<fEventType<<
65c045f0 1441 "Sector="<<isector<<
1442 "Row="<<row<<
1443 "Pad="<<pad<<
1444 "RPad="<<rpad<<
1445 "Max="<<max<<
1446 "Mean="<<mean<<
1447 "RMS="<<rms<<
1448 "GMean="<<gmean<<
1449 "GSigma="<<gsigma<<
03fe1804 1450 "GMeanErr="<<gmeanErr<<
1451 "GSigmaErr="<<gsigmaErr<<
65c045f0 1452 "\n";
1453 if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
1454 }
1455 }
1456}
03fe1804 1457
1458
1459
1460Int_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)
1461{
1462 //
1463 // calculate fourrie transform
1464 // return only frequncies with mag over threshold
1465 // if locMax is spectified only freque with local maxima over theshold is returned
1466
1467 if (! fFFTr2c) return kFALSE;
1468 if (!freq) return kFALSE;
1469
1470 Int_t current=0;
1471 Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
1472 Double_t *in = new Double_t[nPoints];
1473 Double_t *rfft = new Double_t[nPoints];
1474 Double_t *ifft = new Double_t[nPoints];
1475 for (Int_t i=0; i<nPoints; i++){in[i]=input[i];}
1476 fFFTr2c->SetPoints(in);
1477 fFFTr2c->Transform();
1478 fFFTr2c->GetPointsComplex(rfft, ifft);
1479 for (Int_t i=3; i<nPoints/2-3; i++){
1480 Float_t lmag = TMath::Sqrt(rfft[i]*rfft[i]+ifft[i]*ifft[i])/nPoints;
1481 if (lmag<threshold) continue;
1482 if (locMax){
1483 if ( TMath::Sqrt(rfft[i-1]*rfft[i-1]+ifft[i-1]*ifft[i-1])/nPoints>lmag) continue;
1484 if ( TMath::Sqrt(rfft[i+1]*rfft[i+1]+ifft[i+1]*ifft[i+1])/nPoints>lmag) continue;
1485 if ( TMath::Sqrt(rfft[i-2]*rfft[i-2]+ifft[i-2]*ifft[i-2])/nPoints>lmag) continue;
1486 if ( TMath::Sqrt(rfft[i+2]*rfft[i+2]+ifft[i+2]*ifft[i+2])/nPoints>lmag) continue;
1487 if ( TMath::Sqrt(rfft[i-3]*rfft[i-3]+ifft[i-3]*ifft[i-3])/nPoints>lmag) continue;
1488 if ( TMath::Sqrt(rfft[i+3]*rfft[i+3]+ifft[i+3]*ifft[i+3])/nPoints>lmag) continue;
1489 }
1490
1491 freq[current] = Float_t(i)/Float_t(nPoints);
1492 //
1493 re[current] = rfft[i];
1494 im[current] = ifft[i];
1495 mag[current]=lmag;
1496 phi[current]=TMath::ATan2(ifft[i],rfft[i]);
1497 current++;
1498 }
1499 delete [] in;
1500 delete [] rfft;
1501 delete [] ifft;
1502 return current;
1503}
1504