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