]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TPC/AliTPCclustererMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Skip repeating entries (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
... / ...
CommitLineData
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
16/* $Id$ */
17
18//-------------------------------------------------------
19// Implementation of the TPC clusterer
20//
21// 1. The Input data for reconstruction - Options
22// 1.a Simulated data - TTree - invoked Digits2Clusters()
23// 1.b Raw data - Digits2Clusters(AliRawReader* rawReader);
24//
25// 2. The Output data
26// 2.a TTree with clusters - if SetOutput(TTree * tree) invoked
27// 2.b TObjArray - Faster option for HLT
28//
29// 3. Reconstruction setup
30// see AliTPCRecoParam for list of parameters
31// The reconstruction parameterization taken from the
32// AliTPCReconstructor::GetRecoParam()
33// Possible to setup it in reconstruction macro AliTPCReconstructor::SetRecoParam(...)
34//
35//
36//
37// Origin: Marian Ivanov
38//-------------------------------------------------------
39
40#include "Riostream.h"
41#include <TF1.h>
42#include <TFile.h>
43#include <TGraph.h>
44#include <TH1F.h>
45#include <TObjArray.h>
46#include <TRandom.h>
47#include <TTree.h>
48#include <TTreeStream.h>
49
50#include "AliDigits.h"
51#include "AliLoader.h"
52#include "AliLog.h"
53#include "AliMathBase.h"
54#include "AliRawEventHeaderBase.h"
55#include "AliRawReader.h"
56#include "AliRunLoader.h"
57#include "AliSimDigits.h"
58#include "AliTPCCalPad.h"
59#include "AliTPCCalROC.h"
60#include "AliTPCClustersArray.h"
61#include "AliTPCClustersRow.h"
62#include "AliTPCParam.h"
63#include "AliTPCRawStream.h"
64#include "AliTPCRecoParam.h"
65#include "AliTPCReconstructor.h"
66#include "AliTPCcalibDB.h"
67#include "AliTPCclusterInfo.h"
68#include "AliTPCclusterMI.h"
69#include "AliTPCTransform.h"
70#include "AliTPCclustererMI.h"
71
72ClassImp(AliTPCclustererMI)
73
74
75
76AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
77 fBins(0),
78 fSigBins(0),
79 fNSigBins(0),
80 fLoop(0),
81 fMaxBin(0),
82 fMaxTime(0),
83 fMaxPad(0),
84 fSector(-1),
85 fRow(-1),
86 fSign(0),
87 fRx(0),
88 fPadWidth(0),
89 fPadLength(0),
90 fZWidth(0),
91 fPedSubtraction(kFALSE),
92 fEventHeader(0),
93 fTimeStamp(0),
94 fEventType(0),
95 fInput(0),
96 fOutput(0),
97 fOutputArray(0),
98 fRowCl(0),
99 fRowDig(0),
100 fParam(0),
101 fNcluster(0),
102 fDebugStreamer(0),
103 fRecoParam(0),
104 fBDumpSignal(kFALSE)
105{
106 //
107 // COSNTRUCTOR
108 // param - tpc parameters for given file
109 // recoparam - reconstruction parameters
110 //
111 fInput =0;
112 fParam = par;
113 if (recoParam) {
114 fRecoParam = recoParam;
115 }else{
116 //set default parameters if not specified
117 fRecoParam = AliTPCReconstructor::GetRecoParam();
118 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
119 }
120 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
121 Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
122 fRowCl= new AliTPCClustersRow();
123 fRowCl->SetClass("AliTPCclusterMI");
124 fRowCl->SetArray(1);
125
126}
127//______________________________________________________________
128AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
129 :TObject(param),
130 fBins(0),
131 fSigBins(0),
132 fNSigBins(0),
133 fLoop(0),
134 fMaxBin(0),
135 fMaxTime(0),
136 fMaxPad(0),
137 fSector(-1),
138 fRow(-1),
139 fSign(0),
140 fRx(0),
141 fPadWidth(0),
142 fPadLength(0),
143 fZWidth(0),
144 fPedSubtraction(kFALSE),
145 fEventHeader(0),
146 fTimeStamp(0),
147 fEventType(0),
148 fInput(0),
149 fOutput(0),
150 fOutputArray(0),
151 fRowCl(0),
152 fRowDig(0),
153 fParam(0),
154 fNcluster(0),
155 fDebugStreamer(0),
156 fRecoParam(0),
157 fBDumpSignal(kFALSE)
158{
159 //
160 // dummy
161 //
162 fMaxBin = param.fMaxBin;
163}
164//______________________________________________________________
165AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
166{
167 //
168 // assignment operator - dummy
169 //
170 fMaxBin=param.fMaxBin;
171 return (*this);
172}
173//______________________________________________________________
174AliTPCclustererMI::~AliTPCclustererMI(){
175 //
176 //
177 //
178 if (fDebugStreamer) delete fDebugStreamer;
179 if (fOutputArray){
180 //fOutputArray->Delete();
181 delete fOutputArray;
182 }
183}
184
185void AliTPCclustererMI::SetInput(TTree * tree)
186{
187 //
188 // set input tree with digits
189 //
190 fInput = tree;
191 if (!fInput->GetBranch("Segment")){
192 cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
193 fInput=0;
194 return;
195 }
196}
197
198void AliTPCclustererMI::SetOutput(TTree * tree)
199{
200 //
201 // Set the output tree
202 // If not set the ObjArray used - Option for HLT
203 //
204 if (!tree) return;
205 fOutput= tree;
206 AliTPCClustersRow clrow;
207 AliTPCClustersRow *pclrow=&clrow;
208 clrow.SetClass("AliTPCclusterMI");
209 clrow.SetArray(1); // to make Clones array
210 fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
211}
212
213
214void AliTPCclustererMI::FillRow(){
215 //
216 // fill the output container -
217 // 2 Options possible
218 // Tree
219 // TObjArray
220 //
221 if (fOutput) fOutput->Fill();
222 if (!fOutput){
223 //
224 if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
225 if (fRowCl) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
226 }
227}
228
229Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
230 // sigma y2 = in digits - we don't know the angle
231 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
232 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
233 (fPadWidth*fPadWidth);
234 Float_t sres = 0.25;
235 Float_t res = sd2+sres;
236 return res;
237}
238
239
240Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){
241 //sigma z2 = in digits - angle estimated supposing vertex constraint
242 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
243 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
244 Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
245 angular*=angular;
246 angular/=12.;
247 Float_t sres = fParam->GetZSigma()/fZWidth;
248 sres *=sres;
249 Float_t res = angular +sd2+sres;
250 return res;
251}
252
253void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
254AliTPCclusterMI &c)
255{
256 //
257 // k - Make cluster at position k
258 // bins - 2 D array of signals mapped to 1 dimensional array -
259 // max - the number of time bins er one dimension
260 // c - refernce to cluster to be filled
261 //
262 Int_t i0=k/max; //central pad
263 Int_t j0=k%max; //central time bin
264
265 // set pointers to data
266 //Int_t dummy[5] ={0,0,0,0,0};
267 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
268 for (Int_t di=-2;di<=2;di++){
269 matrix[di+2] = &bins[k+di*max];
270 }
271 //build matrix with virtual charge
272 Float_t sigmay2= GetSigmaY2(j0);
273 Float_t sigmaz2= GetSigmaZ2(j0);
274
275 Float_t vmatrix[5][5];
276 vmatrix[2][2] = matrix[2][0];
277 c.SetType(0);
278 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
279 for (Int_t di =-1;di <=1;di++)
280 for (Int_t dj =-1;dj <=1;dj++){
281 Float_t amp = matrix[di+2][dj];
282 if ( (amp<2) && (fLoop<2)){
283 // if under threshold - calculate virtual charge
284 Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
285 amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
286 if (amp>2) amp = 2;
287 vmatrix[2+di][2+dj]=amp;
288 vmatrix[2+2*di][2+2*dj]=0;
289 if ( (di*dj)!=0){
290 //DIAGONAL ELEMENTS
291 vmatrix[2+2*di][2+dj] =0;
292 vmatrix[2+di][2+2*dj] =0;
293 }
294 continue;
295 }
296 if (amp<4){
297 //if small amplitude - below 2 x threshold - don't consider other one
298 vmatrix[2+di][2+dj]=amp;
299 vmatrix[2+2*di][2+2*dj]=0; // don't take to the account next bin
300 if ( (di*dj)!=0){
301 //DIAGONAL ELEMENTS
302 vmatrix[2+2*di][2+dj] =0;
303 vmatrix[2+di][2+2*dj] =0;
304 }
305 continue;
306 }
307 //if bigger then take everything
308 vmatrix[2+di][2+dj]=amp;
309 vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
310 if ( (di*dj)!=0){
311 //DIAGONAL ELEMENTS
312 vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
313 vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
314 }
315 }
316
317
318
319 Float_t sumw=0;
320 Float_t sumiw=0;
321 Float_t sumi2w=0;
322 Float_t sumjw=0;
323 Float_t sumj2w=0;
324 //
325 for (Int_t i=-2;i<=2;i++)
326 for (Int_t j=-2;j<=2;j++){
327 Float_t amp = vmatrix[i+2][j+2];
328
329 sumw += amp;
330 sumiw += i*amp;
331 sumi2w += i*i*amp;
332 sumjw += j*amp;
333 sumj2w += j*j*amp;
334 }
335 //
336 Float_t meani = sumiw/sumw;
337 Float_t mi2 = sumi2w/sumw-meani*meani;
338 Float_t meanj = sumjw/sumw;
339 Float_t mj2 = sumj2w/sumw-meanj*meanj;
340 //
341 Float_t ry = mi2/sigmay2;
342 Float_t rz = mj2/sigmaz2;
343
344 //
345 if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
346 if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) {
347 //
348 //if cluster looks like expected or Unfolding not switched on
349 //standard COG is used
350 //+1.2 deviation from expected sigma accepted
351 // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
352
353 meani +=i0;
354 meanj +=j0;
355 //set cluster parameters
356 c.SetQ(sumw);
357 c.SetPad(meani-2.5);
358 c.SetTimeBin(meanj-3);
359 c.SetSigmaY2(mi2);
360 c.SetSigmaZ2(mj2);
361 c.SetType(0);
362 AddCluster(c,(Float_t*)vmatrix,k);
363 return;
364 }
365 //
366 //unfolding when neccessary
367 //
368
369 Float_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3
370 Float_t dummy[7]={0,0,0,0,0,0};
371 for (Int_t di=-3;di<=3;di++){
372 matrix2[di+3] = &bins[k+di*max];
373 if ((k+di*max)<3) matrix2[di+3] = &dummy[3];
374 if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3];
375 }
376 Float_t vmatrix2[5][5];
377 Float_t sumu;
378 Float_t overlap;
379 UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
380 //
381 // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
382 meani +=i0;
383 meanj +=j0;
384 //set cluster parameters
385 c.SetQ(sumu);
386 c.SetPad(meani-2.5);
387 c.SetTimeBin(meanj-3);
388 c.SetSigmaY2(mi2);
389 c.SetSigmaZ2(mj2);
390 c.SetType(Char_t(overlap)+1);
391 AddCluster(c,(Float_t*)vmatrix,k);
392
393 //unfolding 2
394 meani-=i0;
395 meanj-=j0;
396 if (gDebug>4)
397 printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
398}
399
400
401
402void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
403 Float_t & sumu, Float_t & overlap )
404{
405 //
406 //unfold cluster from input matrix
407 //data corresponding to cluster writen in recmatrix
408 //output meani and meanj
409
410 //take separatelly y and z
411
412 Float_t sum3i[7] = {0,0,0,0,0,0,0};
413 Float_t sum3j[7] = {0,0,0,0,0,0,0};
414
415 for (Int_t k =0;k<7;k++)
416 for (Int_t l = -1; l<=1;l++){
417 sum3i[k]+=matrix2[k][l];
418 sum3j[k]+=matrix2[l+3][k-3];
419 }
420 Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
421 //
422 //unfold y
423 Float_t sum3wi = 0; //charge minus overlap
424 Float_t sum3wio = 0; //full charge
425 Float_t sum3iw = 0; //sum for mean value
426 for (Int_t dk=-1;dk<=1;dk++){
427 sum3wio+=sum3i[dk+3];
428 if (dk==0){
429 sum3wi+=sum3i[dk+3];
430 }
431 else{
432 Float_t ratio =1;
433 if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
434 sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
435 Float_t xm2 = sum3i[-dk+3];
436 Float_t xm1 = sum3i[+3];
437 Float_t x1 = sum3i[2*dk+3];
438 Float_t x2 = sum3i[3*dk+3];
439 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
440 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
441 ratio = w11/(w11+w12);
442 for (Int_t dl=-1;dl<=1;dl++)
443 mratio[dk+1][dl+1] *= ratio;
444 }
445 Float_t amp = sum3i[dk+3]*ratio;
446 sum3wi+=amp;
447 sum3iw+= dk*amp;
448 }
449 }
450 meani = sum3iw/sum3wi;
451 Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
452
453
454
455 //unfold z
456 Float_t sum3wj = 0; //charge minus overlap
457 Float_t sum3wjo = 0; //full charge
458 Float_t sum3jw = 0; //sum for mean value
459 for (Int_t dk=-1;dk<=1;dk++){
460 sum3wjo+=sum3j[dk+3];
461 if (dk==0){
462 sum3wj+=sum3j[dk+3];
463 }
464 else{
465 Float_t ratio =1;
466 if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
467 (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
468 Float_t xm2 = sum3j[-dk+3];
469 Float_t xm1 = sum3j[+3];
470 Float_t x1 = sum3j[2*dk+3];
471 Float_t x2 = sum3j[3*dk+3];
472 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
473 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
474 ratio = w11/(w11+w12);
475 for (Int_t dl=-1;dl<=1;dl++)
476 mratio[dl+1][dk+1] *= ratio;
477 }
478 Float_t amp = sum3j[dk+3]*ratio;
479 sum3wj+=amp;
480 sum3jw+= dk*amp;
481 }
482 }
483 meanj = sum3jw/sum3wj;
484 Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
485 overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
486 sumu = (sum3wj+sum3wi)/2.;
487
488 if (overlap ==3) {
489 //if not overlap detected remove everything
490 for (Int_t di =-2; di<=2;di++)
491 for (Int_t dj =-2; dj<=2;dj++){
492 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
493 }
494 }
495 else{
496 for (Int_t di =-1; di<=1;di++)
497 for (Int_t dj =-1; dj<=1;dj++){
498 Float_t ratio =1;
499 if (mratio[di+1][dj+1]==1){
500 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
501 if (TMath::Abs(di)+TMath::Abs(dj)>1){
502 recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
503 recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
504 }
505 recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
506 }
507 else
508 {
509 //if we have overlap in direction
510 recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
511 if (TMath::Abs(di)+TMath::Abs(dj)>1){
512 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
513 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
514 //
515 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
516 recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
517 }
518 else{
519 ratio = recmatrix[di+2][dj+2]/matrix2[3][0];
520 recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
521 }
522 }
523 }
524 }
525 if (gDebug>4)
526 printf("%f\n", recmatrix[2][2]);
527
528}
529
530Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
531{
532 //
533 // estimate max
534 Float_t sumteor= 0;
535 Float_t sumamp = 0;
536
537 for (Int_t di = -1;di<=1;di++)
538 for (Int_t dj = -1;dj<=1;dj++){
539 if (vmatrix[2+di][2+dj]>2){
540 Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
541 sumteor += teor*vmatrix[2+di][2+dj];
542 sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
543 }
544 }
545 Float_t max = sumamp/sumteor;
546 return max;
547}
548
549void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
550 //
551 //
552 // Transform cluster to the rotated global coordinata
553 // Assign labels to the cluster
554 // add the cluster to the array
555 // for more details - See AliTPCTranform::Transform(x,i,0,1)
556 Float_t meani = c.GetPad();
557 Float_t meanj = c.GetTimeBin();
558
559 Int_t ki = TMath::Nint(meani);
560 if (ki<0) ki=0;
561 if (ki>=fMaxPad) ki = fMaxPad-1;
562 Int_t kj = TMath::Nint(meanj);
563 if (kj<0) kj=0;
564 if (kj>=fMaxTime-3) kj=fMaxTime-4;
565 // ki and kj shifted as integers coordinata
566 if (fRowDig) {
567 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
568 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
569 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
570 }
571
572 c.SetRow(fRow);
573 c.SetDetector(fSector);
574 Float_t s2 = c.GetSigmaY2();
575 Float_t w=fParam->GetPadPitchWidth(fSector);
576 c.SetSigmaY2(s2*w*w);
577 s2 = c.GetSigmaZ2();
578 c.SetSigmaZ2(s2*fZWidth*fZWidth);
579 //
580 //
581 //
582 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
583 if (!transform) {
584 AliFatal("Tranformations not in calibDB");
585 }
586 Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
587 Int_t i[1]={fSector};
588 transform->Transform(x,i,0,1);
589 c.SetX(x[0]);
590 c.SetY(x[1]);
591 c.SetZ(x[2]);
592 //
593 //
594 if (!fRecoParam->GetBYMirror()){
595 if (fSector%36>17){
596 c.SetY(-c.GetY());
597 }
598 }
599
600 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
601 c.SetType(-(c.GetType()+3)); //edge clusters
602 }
603 if (fLoop==2) c.SetType(100);
604
605 TClonesArray * arr = fRowCl->GetArray();
606 AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
607 // if (fRecoParam->DumpSignal() &&matrix ) {
608// Int_t nbins=0;
609// Float_t *graph =0;
610// if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
611// nbins = fMaxTime;
612// graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
613// }
614// AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
615// cl->SetInfo(info);
616// }
617 if (!fRecoParam->DumpSignal()) {
618 cl->SetInfo(0);
619 }
620
621 fNcluster++;
622}
623
624
625//_____________________________________________________________________________
626void AliTPCclustererMI::Digits2Clusters()
627{
628 //-----------------------------------------------------------------
629 // This is a simple cluster finder.
630 //-----------------------------------------------------------------
631
632 if (!fInput) {
633 Error("Digits2Clusters", "input tree not initialised");
634 return;
635 }
636
637 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
638 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
639 AliSimDigits digarr, *dummy=&digarr;
640 fRowDig = dummy;
641 fInput->GetBranch("Segment")->SetAddress(&dummy);
642 Stat_t nentries = fInput->GetEntries();
643
644 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
645
646 Int_t nclusters = 0;
647
648 for (Int_t n=0; n<nentries; n++) {
649 fInput->GetEvent(n);
650 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
651 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
652 continue;
653 }
654 Int_t row = fRow;
655 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
656 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
657 //
658
659 fRowCl->SetID(digarr.GetID());
660 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
661 fRx=fParam->GetPadRowRadii(fSector,row);
662
663
664 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
665 fZWidth = fParam->GetZWidth();
666 if (fSector < kNIS) {
667 fMaxPad = fParam->GetNPadsLow(row);
668 fSign = (fSector < kNIS/2) ? 1 : -1;
669 fPadLength = fParam->GetPadPitchLength(fSector,row);
670 fPadWidth = fParam->GetPadPitchWidth();
671 } else {
672 fMaxPad = fParam->GetNPadsUp(row);
673 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
674 fPadLength = fParam->GetPadPitchLength(fSector,row);
675 fPadWidth = fParam->GetPadPitchWidth();
676 }
677
678
679 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
680 fBins =new Float_t[fMaxBin];
681 fSigBins =new Int_t[fMaxBin];
682 fNSigBins = 0;
683 memset(fBins,0,sizeof(Float_t)*fMaxBin);
684
685 if (digarr.First()) //MI change
686 do {
687 Float_t dig=digarr.CurrentDigit();
688 if (dig<=fParam->GetZeroSup()) continue;
689 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
690 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
691 Int_t bin = i*fMaxTime+j;
692 fBins[bin]=dig/gain;
693 fSigBins[fNSigBins++]=bin;
694 } while (digarr.Next());
695 digarr.ExpandTrackBuffer();
696
697 FindClusters(noiseROC);
698 FillRow();
699 fRowCl->GetArray()->Clear();
700 nclusters+=fNcluster;
701 delete[] fBins;
702 delete[] fSigBins;
703 }
704
705 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
706}
707
708void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
709{
710//-----------------------------------------------------------------
711// This is a cluster finder for the TPC raw data.
712// The method assumes NO ordering of the altro channels.
713// The pedestal subtraction can be switched on and off
714// using an option of the TPC reconstructor
715//-----------------------------------------------------------------
716
717
718 fRowDig = NULL;
719 AliTPCROC * roc = AliTPCROC::Instance();
720 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
721 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
722 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
723 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
724 //
725 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
726 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
727 if (fEventHeader){
728 fTimeStamp = fEventHeader->Get("Timestamp");
729 fEventType = fEventHeader->Get("Type");
730 }
731
732
733 Int_t nclusters = 0;
734
735 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
736 const Int_t kNIS = fParam->GetNInnerSector();
737 const Int_t kNOS = fParam->GetNOuterSector();
738 const Int_t kNS = kNIS + kNOS;
739 fZWidth = fParam->GetZWidth();
740 Int_t zeroSup = fParam->GetZeroSup();
741 //
742 //alocate memory for sector - maximal case
743 //
744 Float_t** allBins = NULL;
745 Int_t** allSigBins = NULL;
746 Int_t* allNSigBins = NULL;
747 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
748 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
749 allBins = new Float_t*[nRowsMax];
750 allSigBins = new Int_t*[nRowsMax];
751 allNSigBins = new Int_t[nRowsMax];
752 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
753 //
754 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
755 allBins[iRow] = new Float_t[maxBin];
756 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
757 allSigBins[iRow] = new Int_t[maxBin];
758 allNSigBins[iRow]=0;
759 }
760 //
761 // Loop over sectors
762 //
763 for(fSector = 0; fSector < kNS; fSector++) {
764
765 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
766 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
767 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
768 //check the presence of the calibration
769 if (!noiseROC ||!pedestalROC ) {
770 AliError(Form("Missing calibration per sector\t%d\n",fSector));
771 continue;
772 }
773 Int_t nRows = 0;
774 Int_t nDDLs = 0, indexDDL = 0;
775 if (fSector < kNIS) {
776 nRows = fParam->GetNRowLow();
777 fSign = (fSector < kNIS/2) ? 1 : -1;
778 nDDLs = 2;
779 indexDDL = fSector * 2;
780 }
781 else {
782 nRows = fParam->GetNRowUp();
783 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
784 nDDLs = 4;
785 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
786 }
787
788 for (Int_t iRow = 0; iRow < nRows; iRow++) {
789 Int_t maxPad;
790 if (fSector < kNIS)
791 maxPad = fParam->GetNPadsLow(iRow);
792 else
793 maxPad = fParam->GetNPadsUp(iRow);
794
795 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
796 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
797 allNSigBins[iRow] = 0;
798 }
799
800 // Loas the raw data for corresponding DDLs
801 rawReader->Reset();
802 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
803 Int_t digCounter=0;
804 // Begin loop over altro data
805 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
806 Float_t gain =1;
807 Int_t lastPad=-1;
808 while (input.Next()) {
809 if (input.GetSector() != fSector)
810 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
811
812
813 Int_t iRow = input.GetRow();
814 if (iRow < 0 || iRow >= nRows){
815 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
816 iRow, 0, nRows -1));
817 continue;
818 }
819 //pad
820 Int_t iPad = input.GetPad();
821 if (iPad < 0 || iPad >= nPadsMax) {
822 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
823 iPad, 0, nPadsMax-1));
824 continue;
825 }
826 if (iPad!=lastPad){
827 gain = gainROC->GetValue(iRow,iPad);
828 lastPad = iPad;
829 }
830 iPad+=3;
831 //time
832 Int_t iTimeBin = input.GetTime();
833 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
834 continue;
835 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
836 iTimeBin, 0, iTimeBin -1));
837 }
838 iTimeBin+=3;
839
840 //signal
841 Float_t signal = input.GetSignal();
842 if (!calcPedestal && signal <= zeroSup) continue;
843 if (!calcPedestal) {
844 Int_t bin = iPad*fMaxTime+iTimeBin;
845 allBins[iRow][bin] = signal/gain;
846 allSigBins[iRow][allNSigBins[iRow]++] = bin;
847 }else{
848 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
849 }
850 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
851
852 // Temporary
853 digCounter++;
854 } // End of the loop over altro data
855 //
856 //
857 //
858 //
859 // Now loop over rows and perform pedestal subtraction
860 if (digCounter==0) continue;
861 // if (calcPedestal) {
862 if (kTRUE) {
863 for (Int_t iRow = 0; iRow < nRows; iRow++) {
864 Int_t maxPad;
865 if (fSector < kNIS)
866 maxPad = fParam->GetNPadsLow(iRow);
867 else
868 maxPad = fParam->GetNPadsUp(iRow);
869
870 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
871 //
872 // Temporary fix for data production - !!!! MARIAN
873 // The noise calibration should take mean and RMS - currently the Gaussian fit used
874 // In case of double peak - the pad should be rejected
875 //
876 // Line mean - if more than given digits over threshold - make a noise calculation
877 // and pedestal substration
878 if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
879 //
880 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
881 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
882 //Float_t pedestal = TMath::Median(fMaxTime, p);
883 Int_t id[3] = {fSector, iRow, iPad-3};
884 // calib values
885 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
886 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
887 Double_t rmsEvent = rmsCalib;
888 Double_t pedestalEvent = pedestalCalib;
889 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
890 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
891 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
892
893 //
894 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
895 Int_t bin = iPad*fMaxTime+iTimeBin;
896 allBins[iRow][bin] -= pedestalEvent;
897 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
898 allBins[iRow][bin] = 0;
899 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
900 allBins[iRow][bin] = 0;
901 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
902 allBins[iRow][bin] = 0;
903 if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
904 allBins[iRow][bin] = 0;
905 if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
906 }
907 }
908 }
909 }
910 // Now loop over rows and find clusters
911 for (fRow = 0; fRow < nRows; fRow++) {
912 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
913 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
914
915 fRx = fParam->GetPadRowRadii(fSector, fRow);
916 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
917 fPadWidth = fParam->GetPadPitchWidth();
918 if (fSector < kNIS)
919 fMaxPad = fParam->GetNPadsLow(fRow);
920 else
921 fMaxPad = fParam->GetNPadsUp(fRow);
922 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
923
924 fBins = allBins[fRow];
925 fSigBins = allSigBins[fRow];
926 fNSigBins = allNSigBins[fRow];
927
928 FindClusters(noiseROC);
929 FillRow();
930 fRowCl->GetArray()->Clear();
931 nclusters += fNcluster;
932 } // End of loop to find clusters
933 } // End of loop over sectors
934
935 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
936 delete [] allBins[iRow];
937 delete [] allSigBins[iRow];
938 }
939 delete [] allBins;
940 delete [] allSigBins;
941 delete [] allNSigBins;
942
943// if (rawReader->GetEventId() && fOutput ){
944// Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
945// }else{
946// Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), nclusters);
947
948// }
949
950}
951
952void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
953{
954
955 //
956 // add virtual charge at the edge
957 //
958 Double_t kMaxDumpSize = 500000;
959 if (!fOutput) {
960 fBDumpSignal =kFALSE;
961 }else{
962 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
963 }
964
965 fNcluster=0;
966 fLoop=1;
967 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
968 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
969 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
970 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
971 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
972 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
973 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
974 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
975 Int_t i = fSigBins[iSig];
976 if (i%fMaxTime<=crtime) continue;
977 Float_t *b = &fBins[i];
978 //absolute custs
979 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
980 //
981 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
982 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
983 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
984 //
985 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
986 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
987 if (!IsMaximum(*b,fMaxTime,b)) continue;
988 //
989 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
990 if (noise>fRecoParam->GetMaxNoise()) continue;
991 // sigma cuts
992 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
993 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
994 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
995
996 AliTPCclusterMI c; // default cosntruction without info
997 Int_t dummy=0;
998 MakeCluster(i, fMaxTime, fBins, dummy,c);
999
1000 //}
1001 }
1002}
1003
1004
1005Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
1006 //
1007 // process signal on given pad - + streaming of additional information in special mode
1008 //
1009 // id[0] - sector
1010 // id[1] - row
1011 // id[2] - pad
1012
1013 //
1014 // ESTIMATE pedestal and the noise
1015 //
1016 const Int_t kPedMax = 100;
1017 Double_t kMaxDebugSize = 5000000.;
1018 Float_t max = 0;
1019 Float_t maxPos = 0;
1020 Int_t median = -1;
1021 Int_t count0 = 0;
1022 Int_t count1 = 0;
1023 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1024 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1025 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
1026 //
1027 UShort_t histo[kPedMax];
1028 memset(histo,0,kPedMax*sizeof(UShort_t));
1029 for (Int_t i=0; i<fMaxTime; i++){
1030 if (signal[i]<=0) continue;
1031 if (signal[i]>max && i>firstBin) {
1032 max = signal[i];
1033 maxPos = i;
1034 }
1035 if (signal[i]>kPedMax-1) continue;
1036 histo[int(signal[i]+0.5)]++;
1037 count0++;
1038 }
1039 //
1040 for (Int_t i=1; i<kPedMax; i++){
1041 if (count1<count0*0.5) median=i;
1042 count1+=histo[i];
1043 }
1044 // truncated mean
1045 //
1046 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1047 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1048 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
1049 //
1050 for (Int_t idelta=1; idelta<10; idelta++){
1051 if (median-idelta<=0) continue;
1052 if (median+idelta>kPedMax) continue;
1053 if (count06<0.6*count1){
1054 count06+=histo[median-idelta];
1055 mean06 +=histo[median-idelta]*(median-idelta);
1056 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1057 count06+=histo[median+idelta];
1058 mean06 +=histo[median+idelta]*(median+idelta);
1059 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1060 }
1061 if (count09<0.9*count1){
1062 count09+=histo[median-idelta];
1063 mean09 +=histo[median-idelta]*(median-idelta);
1064 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1065 count09+=histo[median+idelta];
1066 mean09 +=histo[median+idelta]*(median+idelta);
1067 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1068 }
1069 if (count10<0.95*count1){
1070 count10+=histo[median-idelta];
1071 mean +=histo[median-idelta]*(median-idelta);
1072 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1073 count10+=histo[median+idelta];
1074 mean +=histo[median+idelta]*(median+idelta);
1075 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1076 }
1077 }
1078 if (count10) {
1079 mean /=count10;
1080 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1081 }
1082 if (count06) {
1083 mean06/=count06;
1084 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1085 }
1086 if (count09) {
1087 mean09/=count09;
1088 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1089 }
1090 rmsEvent = rms09;
1091 //
1092 pedestalEvent = median;
1093 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1094 //
1095 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1096 //
1097 // Dump mean signal info
1098 //
1099 (*fDebugStreamer)<<"Signal"<<
1100 "TimeStamp="<<fTimeStamp<<
1101 "EventType="<<fEventType<<
1102 "Sector="<<uid[0]<<
1103 "Row="<<uid[1]<<
1104 "Pad="<<uid[2]<<
1105 "Max="<<max<<
1106 "MaxPos="<<maxPos<<
1107 //
1108 "Median="<<median<<
1109 "Mean="<<mean<<
1110 "RMS="<<rms<<
1111 "Mean06="<<mean06<<
1112 "RMS06="<<rms06<<
1113 "Mean09="<<mean09<<
1114 "RMS09="<<rms09<<
1115 "RMSCalib="<<rmsCalib<<
1116 "PedCalib="<<pedestalCalib<<
1117 "\n";
1118 //
1119 // fill pedestal histogram
1120 //
1121 AliTPCROC * roc = AliTPCROC::Instance();
1122
1123 //
1124 //
1125 //
1126 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1127 Float_t *dsignal = new Float_t[nchannels];
1128 Float_t *dtime = new Float_t[nchannels];
1129 for (Int_t i=0; i<nchannels; i++){
1130 dtime[i] = i;
1131 dsignal[i] = signal[i];
1132 }
1133
1134 TGraph * graph;
1135 //
1136 // Big signals dumping
1137 //
1138 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1139 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1140 "TimeStamp="<<fTimeStamp<<
1141 "EventType="<<fEventType<<
1142 "Sector="<<uid[0]<<
1143 "Row="<<uid[1]<<
1144 "Pad="<<uid[2]<<
1145 "Graph="<<graph<<
1146 "Max="<<max<<
1147 "MaxPos="<<maxPos<<
1148 //
1149 "Median="<<median<<
1150 "Mean="<<mean<<
1151 "RMS="<<rms<<
1152 "Mean06="<<mean06<<
1153 "RMS06="<<rms06<<
1154 "Mean09="<<mean09<<
1155 "RMS09="<<rms09<<
1156 "\n";
1157 delete graph;
1158
1159 delete [] dsignal;
1160 delete [] dtime;
1161 if (rms06>fRecoParam->GetMaxNoise()) {
1162 pedestalEvent+=1024.;
1163 return 1024+median; // sign noisy channel in debug mode
1164 }
1165 return median;
1166}
1167
1168
1169
1170
1171