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