]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCclustererMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding stream level and debug streamer to the base class
[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//
ebe95b38 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//
1c53abe2 37// Origin: Marian Ivanov
38//-------------------------------------------------------
39
a1e17193 40#include "Riostream.h"
41#include <TF1.h>
1c53abe2 42#include <TFile.h>
a1e17193 43#include <TGraph.h>
44#include <TH1F.h>
45#include <TObjArray.h>
46#include <TRandom.h>
47#include <TTree.h>
48#include <TTreeStream.h>
65c045f0 49
1c53abe2 50#include "AliDigits.h"
a1e17193 51#include "AliLoader.h"
52#include "AliLog.h"
53#include "AliMathBase.h"
ef4ad662 54#include "AliRawEventHeaderBase.h"
a1e17193 55#include "AliRawReader.h"
f8aae377 56#include "AliRunLoader.h"
a1e17193 57#include "AliSimDigits.h"
13116aec 58#include "AliTPCCalPad.h"
59#include "AliTPCCalROC.h"
a1e17193 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"
022ee144 69#include "AliTPCTransform.h"
a1e17193 70#include "AliTPCclustererMI.h"
13116aec 71
1c53abe2 72ClassImp(AliTPCclustererMI)
73
74
75
2f32844a 76AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
77 fBins(0),
5b33a7f5 78 fSigBins(0),
79 fNSigBins(0),
2f32844a 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),
ef4ad662 93 fEventHeader(0),
94 fTimeStamp(0),
95 fEventType(0),
2f32844a 96 fInput(0),
97 fOutput(0),
ebe95b38 98 fOutputArray(0),
2f32844a 99 fRowCl(0),
100 fRowDig(0),
101 fParam(0),
102 fNcluster(0),
2f32844a 103 fDebugStreamer(0),
03fe1804 104 fRecoParam(0),
a525cc34 105 fBDumpSignal(kFALSE)
1c53abe2 106{
97f127bb 107 //
108 // COSNTRUCTOR
109 // param - tpc parameters for given file
110 // recoparam - reconstruction parameters
111 //
22c352f8 112 fIsOldRCUFormat = kFALSE;
1c53abe2 113 fInput =0;
f8aae377 114 fParam = par;
194b0609 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 }
65c045f0 122 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
03fe1804 123 Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
65c045f0 124}
e046d791 125//______________________________________________________________
126AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
127 :TObject(param),
128 fBins(0),
5b33a7f5 129 fSigBins(0),
130 fNSigBins(0),
e046d791 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),
ebe95b38 149 fOutputArray(0),
e046d791 150 fRowCl(0),
151 fRowDig(0),
152 fParam(0),
153 fNcluster(0),
e046d791 154 fDebugStreamer(0),
5b33a7f5 155 fRecoParam(0),
a525cc34 156 fBDumpSignal(kFALSE)
e046d791 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//______________________________________________________________
65c045f0 173AliTPCclustererMI::~AliTPCclustererMI(){
ebe95b38 174 //
175 //
176 //
65c045f0 177 if (fDebugStreamer) delete fDebugStreamer;
ebe95b38 178 if (fOutputArray){
179 fOutputArray->Delete();
180 delete fOutputArray;
181 }
1c53abe2 182}
22c352f8 183
1c53abe2 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 //
ebe95b38 200 // Set the output tree
201 // If not set the ObjArray used - Option for HLT
1c53abe2 202 //
ebe95b38 203 if (!tree) return;
204 fOutput= tree;
1c53abe2 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
ebe95b38 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
1c53abe2 228Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
229 // sigma y2 = in digits - we don't know the angle
753797ce 230 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 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
753797ce 241 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 242 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
a1ec4d07 243 Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
1c53abe2 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
13116aec 252void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
1c53abe2 253AliTPCclusterMI &c)
254{
97f127bb 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 //
1c53abe2 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};
13116aec 266 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
1c53abe2 267 for (Int_t di=-2;di<=2;di++){
268 matrix[di+2] = &bins[k+di*max];
1c53abe2 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);
6c024a0e 277 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
1c53abe2 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;
194b0609 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
1c53abe2 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);
022ee144 356 c.SetPad(meani-2.5);
5bd7ed29 357 c.SetTimeBin(meanj-3);
1c53abe2 358 c.SetSigmaY2(mi2);
359 c.SetSigmaZ2(mj2);
022ee144 360 c.SetType(0);
03fe1804 361 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 362 return;
363 }
364 //
365 //unfolding when neccessary
366 //
367
13116aec 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};
1c53abe2 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);
022ee144 385 c.SetPad(meani-2.5);
386 c.SetTimeBin(meanj-3);
1c53abe2 387 c.SetSigmaY2(mi2);
388 c.SetSigmaZ2(mj2);
389 c.SetType(Char_t(overlap)+1);
03fe1804 390 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 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
13116aec 401void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
1c53abe2 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];
cc5e9db0 438 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
439 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 440 ratio = w11/(w11+w12);
441 for (Int_t dl=-1;dl<=1;dl++)
442 mratio[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];
cc5e9db0 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.);
1c53abe2 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){
cc5e9db0 511 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
1c53abe2 512 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
513 //
cc5e9db0 514 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
1c53abe2 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
03fe1804 548void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
1c53abe2 549 //
1c53abe2 550 //
022ee144 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();
1c53abe2 557
022ee144 558 Int_t ki = TMath::Nint(meani);
1c53abe2 559 if (ki<0) ki=0;
560 if (ki>=fMaxPad) ki = fMaxPad-1;
022ee144 561 Int_t kj = TMath::Nint(meanj);
1c53abe2 562 if (kj<0) kj=0;
563 if (kj>=fMaxTime-3) kj=fMaxTime-4;
022ee144 564 // ki and kj shifted as integers coordinata
f8aae377 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 }
1c53abe2 570
022ee144 571 c.SetRow(fRow);
572 c.SetDetector(fSector);
1c53abe2 573 Float_t s2 = c.GetSigmaY2();
574 Float_t w=fParam->GetPadPitchWidth(fSector);
1c53abe2 575 c.SetSigmaY2(s2*w*w);
576 s2 = c.GetSigmaZ2();
022ee144 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 //
adefcaa6 593 if (!fRecoParam->GetBYMirror()){
594 if (fSector%36>17){
022ee144 595 c.SetY(-c.GetY());
adefcaa6 596 }
597 }
508541c7 598
1c53abe2 599 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
1c53abe2 600 c.SetType(-(c.GetType()+3)); //edge clusters
601 }
602 if (fLoop==2) c.SetType(100);
603
604 TClonesArray * arr = fRowCl->GetArray();
03fe1804 605 AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
b127a65f 606 if (fRecoParam->DumpSignal() &&matrix ) {
03fe1804 607 Int_t nbins=0;
608 Float_t *graph =0;
7041b196 609 if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
03fe1804 610 nbins = fMaxTime;
611 graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
612 }
613 AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
614 cl->SetInfo(info);
615 }
b127a65f 616 if (!fRecoParam->DumpSignal()) {
617 cl->SetInfo(0);
618 }
1c53abe2 619
620 fNcluster++;
621}
622
623
624//_____________________________________________________________________________
f8aae377 625void AliTPCclustererMI::Digits2Clusters()
1c53abe2 626{
627 //-----------------------------------------------------------------
628 // This is a simple cluster finder.
629 //-----------------------------------------------------------------
1c53abe2 630
f8aae377 631 if (!fInput) {
632 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 633 return;
634 }
635
13116aec 636 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 637 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1c53abe2 638 AliSimDigits digarr, *dummy=&digarr;
639 fRowDig = dummy;
640 fInput->GetBranch("Segment")->SetAddress(&dummy);
641 Stat_t nentries = fInput->GetEntries();
642
daac2a58 643 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 644
1c53abe2 645 Int_t nclusters = 0;
13116aec 646
1c53abe2 647 for (Int_t n=0; n<nentries; n++) {
648 fInput->GetEvent(n);
508541c7 649 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 650 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
651 continue;
652 }
508541c7 653 Int_t row = fRow;
13116aec 654 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 655 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 656 //
ebe95b38 657 fRowCl= new AliTPCClustersRow();
658 fRowCl->SetClass("AliTPCclusterMI");
659 fRowCl->SetArray(1);
1c53abe2 660
ebe95b38 661 fRowCl->SetID(digarr.GetID());
662 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
f8aae377 663 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 664
665
f8aae377 666 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 667 fZWidth = fParam->GetZWidth();
668 if (fSector < kNIS) {
f8aae377 669 fMaxPad = fParam->GetNPadsLow(row);
1c53abe2 670 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 671 fPadLength = fParam->GetPadPitchLength(fSector,row);
672 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 673 } else {
f8aae377 674 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 675 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 676 fPadLength = fParam->GetPadPitchLength(fSector,row);
677 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 678 }
679
680
681 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 682 fBins =new Float_t[fMaxBin];
5b33a7f5 683 fSigBins =new Int_t[fMaxBin];
684 fNSigBins = 0;
13116aec 685 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 686
687 if (digarr.First()) //MI change
688 do {
13116aec 689 Float_t dig=digarr.CurrentDigit();
f8aae377 690 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 691 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 692 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 693 Int_t bin = i*fMaxTime+j;
694 fBins[bin]=dig/gain;
695 fSigBins[fNSigBins++]=bin;
1c53abe2 696 } while (digarr.Next());
697 digarr.ExpandTrackBuffer();
698
940ed1f0 699 FindClusters(noiseROC);
ebe95b38 700 FillRow();
701 delete fRowCl;
88cb7938 702 nclusters+=fNcluster;
5b33a7f5 703 delete[] fBins;
704 delete[] fSigBins;
88cb7938 705 }
f8aae377 706
19dd5b2f 707 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
f8aae377 708}
709
710void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
711{
712//-----------------------------------------------------------------
22c352f8 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
f8aae377 717//-----------------------------------------------------------------
718
f8aae377 719
f8aae377 720 fRowDig = NULL;
adefcaa6 721 AliTPCROC * roc = AliTPCROC::Instance();
22c352f8 722 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 723 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
724 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
d6834f5f 725 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
726 //
727 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 728 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
729 if (fEventHeader){
730 fTimeStamp = fEventHeader->Get("Timestamp");
731 fEventType = fEventHeader->Get("Type");
732 }
733
22c352f8 734
f8aae377 735 Int_t nclusters = 0;
88cb7938 736
daac2a58 737 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 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();
adefcaa6 743 //
744 //alocate memory for sector - maximal case
745 //
22c352f8 746 Float_t** allBins = NULL;
5b33a7f5 747 Int_t** allSigBins = NULL;
748 Int_t* allNSigBins = NULL;
adefcaa6 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];
5b33a7f5 752 allSigBins = new Int_t*[nRowsMax];
753 allNSigBins = new Int_t[nRowsMax];
adefcaa6 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];
adefcaa6 758 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 759 allSigBins[iRow] = new Int_t[maxBin];
760 allNSigBins[iRow]=0;
adefcaa6 761 }
762 //
22c352f8 763 // Loop over sectors
adefcaa6 764 //
22c352f8 765 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 766
940ed1f0 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
b24bd339 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 }
22c352f8 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
22c352f8 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
22c352f8 798 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 799 allNSigBins[iRow] = 0;
22c352f8 800 }
801
802 // Loas the raw data for corresponding DDLs
803 rawReader->Reset();
22c352f8 804 input.SetOldRCUFormat(fIsOldRCUFormat);
362c9d61 805 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
adefcaa6 806 Int_t digCounter=0;
22c352f8 807 // Begin loop over altro data
adefcaa6 808 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
809 Float_t gain =1;
810 Int_t lastPad=-1;
22c352f8 811 while (input.Next()) {
22c352f8 812 if (input.GetSector() != fSector)
813 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
814
22c352f8 815
816 Int_t iRow = input.GetRow();
3f0af4ba 817 if (iRow < 0 || iRow >= nRows){
818 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 819 iRow, 0, nRows -1));
3f0af4ba 820 continue;
821 }
adefcaa6 822 //pad
823 Int_t iPad = input.GetPad();
3f0af4ba 824 if (iPad < 0 || iPad >= nPadsMax) {
825 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 826 iPad, 0, nPadsMax-1));
3f0af4ba 827 continue;
828 }
adefcaa6 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();
daac2a58 836 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
837 continue;
22c352f8 838 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 839 iTimeBin, 0, iTimeBin -1));
daac2a58 840 }
adefcaa6 841 iTimeBin+=3;
3f0af4ba 842
adefcaa6 843 //signal
22c352f8 844 Float_t signal = input.GetSignal();
adefcaa6 845 if (!calcPedestal && signal <= zeroSup) continue;
940ed1f0 846 if (!calcPedestal) {
5b33a7f5 847 Int_t bin = iPad*fMaxTime+iTimeBin;
848 allBins[iRow][bin] = signal/gain;
849 allSigBins[iRow][allNSigBins[iRow]++] = bin;
940ed1f0 850 }else{
851 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
852 }
aba7fdfc 853 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 854
855 // Temporary
856 digCounter++;
22c352f8 857 } // End of the loop over altro data
adefcaa6 858 //
859 //
aba7fdfc 860 //
861 //
22c352f8 862 // Now loop over rows and perform pedestal subtraction
adefcaa6 863 if (digCounter==0) continue;
aba7fdfc 864 // if (calcPedestal) {
865 if (kTRUE) {
22c352f8 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
940ed1f0 873 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
aba7fdfc 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 //
940ed1f0 883 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 884 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 885 //Float_t pedestal = TMath::Median(fMaxTime, p);
886 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 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 //
22c352f8 897 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
5b33a7f5 898 Int_t bin = iPad*fMaxTime+iTimeBin;
899 allBins[iRow][bin] -= pedestalEvent;
7fef31a6 900 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
5b33a7f5 901 allBins[iRow][bin] = 0;
7fef31a6 902 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
5b33a7f5 903 allBins[iRow][bin] = 0;
22c352f8 904 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
5b33a7f5 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;
22c352f8 909 }
910 }
f8aae377 911 }
22c352f8 912 }
22c352f8 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));
ebe95b38 919 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
22c352f8 920
921 fRx = fParam->GetPadRowRadii(fSector, fRow);
922 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 923 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 924 if (fSector < kNIS)
925 fMaxPad = fParam->GetNPadsLow(fRow);
926 else
927 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 928 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 929
930 fBins = allBins[fRow];
5b33a7f5 931 fSigBins = allSigBins[fRow];
932 fNSigBins = allNSigBins[fRow];
22c352f8 933
940ed1f0 934 FindClusters(noiseROC);
ebe95b38 935 FillRow();
22c352f8 936 delete fRowCl;
937 nclusters += fNcluster;
938 } // End of loop to find clusters
22c352f8 939 } // End of loop over sectors
adefcaa6 940
941 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
942 delete [] allBins[iRow];
5b33a7f5 943 delete [] allSigBins[iRow];
adefcaa6 944 }
945 delete [] allBins;
5b33a7f5 946 delete [] allSigBins;
947 delete [] allNSigBins;
adefcaa6 948
d09996c4 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);
ebe95b38 953
d09996c4 954// }
ebe95b38 955
f8aae377 956}
957
940ed1f0 958void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 959{
940ed1f0 960
961 //
962 // add virtual charge at the edge
963 //
7041b196 964 Double_t kMaxDumpSize = 500000;
ebe95b38 965 if (!fOutput) {
966 fBDumpSignal =kFALSE;
967 }else{
968 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
969 }
970
f8aae377 971 fNcluster=0;
f8aae377 972 fLoop=1;
a1ec4d07 973 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 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();
5b33a7f5 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];
940ed1f0 984 //absolute custs
03fe1804 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
aba7fdfc 988 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
989 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 990 //
991 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 992 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 993 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 994 //
995 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 996 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 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
03fe1804 1002 AliTPCclusterMI c(kFALSE); // default cosntruction without info
f8aae377 1003 Int_t dummy=0;
1004 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1005
f8aae377 1006 //}
1007 }
8569a2b0 1008}
65c045f0 1009
1010
940ed1f0 1011Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1012 //
1013 // process signal on given pad - + streaming of additional information in special mode
1014 //
1015 // id[0] - sector
1016 // id[1] - row
adefcaa6 1017 // id[2] - pad
1018
1019 //
1020 // ESTIMATE pedestal and the noise
1021 //
adefcaa6 1022 const Int_t kPedMax = 100;
7041b196 1023 Double_t kMaxDebugSize = 5000000.;
adefcaa6 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;
940ed1f0 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();
65c045f0 1032 //
adefcaa6 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;
940ed1f0 1037 if (signal[i]>max && i>firstBin) {
65c045f0 1038 max = signal[i];
1039 maxPos = i;
adefcaa6 1040 }
1041 if (signal[i]>kPedMax-1) continue;
1042 histo[int(signal[i]+0.5)]++;
1043 count0++;
65c045f0 1044 }
7fef31a6 1045 //
adefcaa6 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
65c045f0 1051 //
7041b196 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 ;
adefcaa6 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);
65c045f0 1074 }
adefcaa6 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 }
3f0af4ba 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 }
940ed1f0 1096 rmsEvent = rms09;
adefcaa6 1097 //
940ed1f0 1098 pedestalEvent = median;
adefcaa6 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"<<
ef4ad662 1106 "TimeStamp="<<fTimeStamp<<
1107 "EventType="<<fEventType<<
adefcaa6 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<<
940ed1f0 1121 "RMSCalib="<<rmsCalib<<
1122 "PedCalib="<<pedestalCalib<<
adefcaa6 1123 "\n";
1124 //
1125 // fill pedestal histogram
1126 //
1127 AliTPCROC * roc = AliTPCROC::Instance();
a525cc34 1128
65c045f0 1129 //
adefcaa6 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 }
03fe1804 1139
65c045f0 1140 TGraph * graph;
7fef31a6 1141 //
a525cc34 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;
7fef31a6 1164
65c045f0 1165 delete [] dsignal;
1166 delete [] dtime;
940ed1f0 1167 if (rms06>fRecoParam->GetMaxNoise()) {
1168 pedestalEvent+=1024.;
1169 return 1024+median; // sign noisy channel in debug mode
1170 }
65c045f0 1171 return median;
1172}
1173
1174
1175
03fe1804 1176
03fe1804 1177