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