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