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