]> 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
- Small fix for CAF that messed-up user tasks after filters
[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),
ef4ad662 92 fEventHeader(0),
93 fTimeStamp(0),
94 fEventType(0),
2f32844a 95 fInput(0),
96 fOutput(0),
ebe95b38 97 fOutputArray(0),
2f32844a 98 fRowCl(0),
99 fRowDig(0),
100 fParam(0),
101 fNcluster(0),
2f32844a 102 fDebugStreamer(0),
03fe1804 103 fRecoParam(0),
a525cc34 104 fBDumpSignal(kFALSE)
1c53abe2 105{
97f127bb 106 //
107 // COSNTRUCTOR
108 // param - tpc parameters for given file
109 // recoparam - reconstruction parameters
110 //
1c53abe2 111 fInput =0;
f8aae377 112 fParam = par;
194b0609 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 }
65c045f0 120 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
03fe1804 121 Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
d101caf3 122 fRowCl= new AliTPCClustersRow();
123 fRowCl->SetClass("AliTPCclusterMI");
124 fRowCl->SetArray(1);
125
65c045f0 126}
e046d791 127//______________________________________________________________
128AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
129 :TObject(param),
130 fBins(0),
5b33a7f5 131 fSigBins(0),
132 fNSigBins(0),
e046d791 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),
e046d791 145 fEventHeader(0),
146 fTimeStamp(0),
147 fEventType(0),
148 fInput(0),
149 fOutput(0),
ebe95b38 150 fOutputArray(0),
e046d791 151 fRowCl(0),
152 fRowDig(0),
153 fParam(0),
154 fNcluster(0),
e046d791 155 fDebugStreamer(0),
5b33a7f5 156 fRecoParam(0),
a525cc34 157 fBDumpSignal(kFALSE)
e046d791 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//______________________________________________________________
65c045f0 174AliTPCclustererMI::~AliTPCclustererMI(){
ebe95b38 175 //
176 //
177 //
65c045f0 178 if (fDebugStreamer) delete fDebugStreamer;
ebe95b38 179 if (fOutputArray){
af86c1fd 180 //fOutputArray->Delete();
ebe95b38 181 delete fOutputArray;
182 }
1c53abe2 183}
22c352f8 184
1c53abe2 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 //
ebe95b38 201 // Set the output tree
202 // If not set the ObjArray used - Option for HLT
1c53abe2 203 //
ebe95b38 204 if (!tree) return;
205 fOutput= tree;
1c53abe2 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
ebe95b38 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 //
af86c1fd 224 if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
ebe95b38 225 if (fRowCl) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
226 }
227}
228
1c53abe2 229Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
230 // sigma y2 = in digits - we don't know the angle
753797ce 231 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 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
753797ce 242 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
1c53abe2 243 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
a1ec4d07 244 Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
1c53abe2 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
13116aec 253void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
1c53abe2 254AliTPCclusterMI &c)
255{
97f127bb 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 //
1c53abe2 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};
13116aec 267 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
1c53abe2 268 for (Int_t di=-2;di<=2;di++){
269 matrix[di+2] = &bins[k+di*max];
1c53abe2 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);
6c024a0e 278 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
1c53abe2 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;
194b0609 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
1c53abe2 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);
022ee144 357 c.SetPad(meani-2.5);
5bd7ed29 358 c.SetTimeBin(meanj-3);
1c53abe2 359 c.SetSigmaY2(mi2);
360 c.SetSigmaZ2(mj2);
022ee144 361 c.SetType(0);
03fe1804 362 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 363 return;
364 }
365 //
366 //unfolding when neccessary
367 //
368
13116aec 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};
1c53abe2 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);
022ee144 386 c.SetPad(meani-2.5);
387 c.SetTimeBin(meanj-3);
1c53abe2 388 c.SetSigmaY2(mi2);
389 c.SetSigmaZ2(mj2);
390 c.SetType(Char_t(overlap)+1);
03fe1804 391 AddCluster(c,(Float_t*)vmatrix,k);
1c53abe2 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
13116aec 402void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
1c53abe2 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];
cc5e9db0 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.);
1c53abe2 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];
cc5e9db0 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.);
1c53abe2 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){
cc5e9db0 512 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
1c53abe2 513 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
514 //
cc5e9db0 515 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
1c53abe2 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
03fe1804 549void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
1c53abe2 550 //
1c53abe2 551 //
022ee144 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();
1c53abe2 558
022ee144 559 Int_t ki = TMath::Nint(meani);
1c53abe2 560 if (ki<0) ki=0;
561 if (ki>=fMaxPad) ki = fMaxPad-1;
022ee144 562 Int_t kj = TMath::Nint(meanj);
1c53abe2 563 if (kj<0) kj=0;
564 if (kj>=fMaxTime-3) kj=fMaxTime-4;
022ee144 565 // ki and kj shifted as integers coordinata
f8aae377 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 }
1c53abe2 571
022ee144 572 c.SetRow(fRow);
573 c.SetDetector(fSector);
1c53abe2 574 Float_t s2 = c.GetSigmaY2();
575 Float_t w=fParam->GetPadPitchWidth(fSector);
1c53abe2 576 c.SetSigmaY2(s2*w*w);
577 s2 = c.GetSigmaZ2();
022ee144 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 //
adefcaa6 594 if (!fRecoParam->GetBYMirror()){
595 if (fSector%36>17){
022ee144 596 c.SetY(-c.GetY());
adefcaa6 597 }
598 }
508541c7 599
1c53abe2 600 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
1c53abe2 601 c.SetType(-(c.GetType()+3)); //edge clusters
602 }
603 if (fLoop==2) c.SetType(100);
604
605 TClonesArray * arr = fRowCl->GetArray();
03fe1804 606 AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
d101caf3 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// }
b127a65f 617 if (!fRecoParam->DumpSignal()) {
618 cl->SetInfo(0);
619 }
1c53abe2 620
621 fNcluster++;
622}
623
624
625//_____________________________________________________________________________
f8aae377 626void AliTPCclustererMI::Digits2Clusters()
1c53abe2 627{
628 //-----------------------------------------------------------------
629 // This is a simple cluster finder.
630 //-----------------------------------------------------------------
1c53abe2 631
f8aae377 632 if (!fInput) {
633 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 634 return;
635 }
636
13116aec 637 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 638 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1c53abe2 639 AliSimDigits digarr, *dummy=&digarr;
640 fRowDig = dummy;
641 fInput->GetBranch("Segment")->SetAddress(&dummy);
642 Stat_t nentries = fInput->GetEntries();
643
daac2a58 644 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 645
1c53abe2 646 Int_t nclusters = 0;
13116aec 647
1c53abe2 648 for (Int_t n=0; n<nentries; n++) {
649 fInput->GetEvent(n);
508541c7 650 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 651 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
652 continue;
653 }
508541c7 654 Int_t row = fRow;
13116aec 655 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 656 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 657 //
1c53abe2 658
ebe95b38 659 fRowCl->SetID(digarr.GetID());
660 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
f8aae377 661 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 662
663
f8aae377 664 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 665 fZWidth = fParam->GetZWidth();
666 if (fSector < kNIS) {
f8aae377 667 fMaxPad = fParam->GetNPadsLow(row);
1c53abe2 668 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 669 fPadLength = fParam->GetPadPitchLength(fSector,row);
670 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 671 } else {
f8aae377 672 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 673 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 674 fPadLength = fParam->GetPadPitchLength(fSector,row);
675 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 676 }
677
678
679 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 680 fBins =new Float_t[fMaxBin];
5b33a7f5 681 fSigBins =new Int_t[fMaxBin];
682 fNSigBins = 0;
13116aec 683 memset(fBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 684
685 if (digarr.First()) //MI change
686 do {
13116aec 687 Float_t dig=digarr.CurrentDigit();
f8aae377 688 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 689 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 690 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
5b33a7f5 691 Int_t bin = i*fMaxTime+j;
692 fBins[bin]=dig/gain;
693 fSigBins[fNSigBins++]=bin;
1c53abe2 694 } while (digarr.Next());
695 digarr.ExpandTrackBuffer();
696
940ed1f0 697 FindClusters(noiseROC);
ebe95b38 698 FillRow();
d101caf3 699 fRowCl->GetArray()->Clear();
88cb7938 700 nclusters+=fNcluster;
5b33a7f5 701 delete[] fBins;
702 delete[] fSigBins;
88cb7938 703 }
f8aae377 704
19dd5b2f 705 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
f8aae377 706}
707
708void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
709{
710//-----------------------------------------------------------------
22c352f8 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
f8aae377 715//-----------------------------------------------------------------
716
f8aae377 717
f8aae377 718 fRowDig = NULL;
adefcaa6 719 AliTPCROC * roc = AliTPCROC::Instance();
22c352f8 720 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 721 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
722 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
d6834f5f 723 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
724 //
725 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
ef4ad662 726 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
727 if (fEventHeader){
728 fTimeStamp = fEventHeader->Get("Timestamp");
729 fEventType = fEventHeader->Get("Type");
730 }
731
22c352f8 732
f8aae377 733 Int_t nclusters = 0;
88cb7938 734
daac2a58 735 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
f8aae377 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();
adefcaa6 741 //
742 //alocate memory for sector - maximal case
743 //
22c352f8 744 Float_t** allBins = NULL;
5b33a7f5 745 Int_t** allSigBins = NULL;
746 Int_t* allNSigBins = NULL;
adefcaa6 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];
5b33a7f5 750 allSigBins = new Int_t*[nRowsMax];
751 allNSigBins = new Int_t[nRowsMax];
adefcaa6 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];
adefcaa6 756 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 757 allSigBins[iRow] = new Int_t[maxBin];
758 allNSigBins[iRow]=0;
adefcaa6 759 }
760 //
22c352f8 761 // Loop over sectors
adefcaa6 762 //
22c352f8 763 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 764
940ed1f0 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
b24bd339 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 }
22c352f8 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
22c352f8 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
22c352f8 796 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
5b33a7f5 797 allNSigBins[iRow] = 0;
22c352f8 798 }
799
800 // Loas the raw data for corresponding DDLs
801 rawReader->Reset();
362c9d61 802 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
adefcaa6 803 Int_t digCounter=0;
22c352f8 804 // Begin loop over altro data
adefcaa6 805 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
806 Float_t gain =1;
807 Int_t lastPad=-1;
22c352f8 808 while (input.Next()) {
22c352f8 809 if (input.GetSector() != fSector)
810 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
811
22c352f8 812
813 Int_t iRow = input.GetRow();
9546a7ff 814 if (iRow < 0){
815 continue;
816 }
817
3f0af4ba 818 if (iRow < 0 || iRow >= nRows){
819 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
22c352f8 820 iRow, 0, nRows -1));
3f0af4ba 821 continue;
822 }
adefcaa6 823 //pad
824 Int_t iPad = input.GetPad();
3f0af4ba 825 if (iPad < 0 || iPad >= nPadsMax) {
826 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 827 iPad, 0, nPadsMax-1));
3f0af4ba 828 continue;
829 }
adefcaa6 830 if (iPad!=lastPad){
831 gain = gainROC->GetValue(iRow,iPad);
832 lastPad = iPad;
833 }
834 iPad+=3;
835 //time
836 Int_t iTimeBin = input.GetTime();
daac2a58 837 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
838 continue;
22c352f8 839 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 840 iTimeBin, 0, iTimeBin -1));
daac2a58 841 }
adefcaa6 842 iTimeBin+=3;
3f0af4ba 843
adefcaa6 844 //signal
22c352f8 845 Float_t signal = input.GetSignal();
adefcaa6 846 if (!calcPedestal && signal <= zeroSup) continue;
9546a7ff 847 if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
848 Double_t x[]={iRow,iPad,iTimeBin};
849 Int_t i[]={fSector};
850 AliTPCTransform trafo;
851 trafo.Transform(x,i,0,1);
852 Double_t gx[3]={x[0],x[1],x[2]};
853 trafo.RotatedGlobal2Global(fSector,gx);
854
855 (*fDebugStreamer)<<"Digits"<<
856 "sec="<<fSector<<
857 "row="<<iRow<<
858 "pad="<<iPad<<
859 "sig="<<signal<<
860 "x="<<x[0]<<
861 "y="<<x[1]<<
862 "z="<<x[2]<<
863 "gx="<<gx[0]<<
864 "gy="<<gx[1]<<
865 "gz="<<gx[2]<<
866 "\n";
867 }
868
869
940ed1f0 870 if (!calcPedestal) {
5b33a7f5 871 Int_t bin = iPad*fMaxTime+iTimeBin;
872 allBins[iRow][bin] = signal/gain;
873 allSigBins[iRow][allNSigBins[iRow]++] = bin;
940ed1f0 874 }else{
875 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
876 }
aba7fdfc 877 allBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
3f0af4ba 878
879 // Temporary
880 digCounter++;
22c352f8 881 } // End of the loop over altro data
adefcaa6 882 //
883 //
aba7fdfc 884 //
885 //
22c352f8 886 // Now loop over rows and perform pedestal subtraction
adefcaa6 887 if (digCounter==0) continue;
aba7fdfc 888 // if (calcPedestal) {
889 if (kTRUE) {
22c352f8 890 for (Int_t iRow = 0; iRow < nRows; iRow++) {
891 Int_t maxPad;
892 if (fSector < kNIS)
893 maxPad = fParam->GetNPadsLow(iRow);
894 else
895 maxPad = fParam->GetNPadsUp(iRow);
896
940ed1f0 897 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
aba7fdfc 898 //
899 // Temporary fix for data production - !!!! MARIAN
900 // The noise calibration should take mean and RMS - currently the Gaussian fit used
901 // In case of double peak - the pad should be rejected
902 //
903 // Line mean - if more than given digits over threshold - make a noise calculation
904 // and pedestal substration
905 if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
906 //
940ed1f0 907 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 908 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 909 //Float_t pedestal = TMath::Median(fMaxTime, p);
910 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 911 // calib values
912 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
913 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
914 Double_t rmsEvent = rmsCalib;
915 Double_t pedestalEvent = pedestalCalib;
916 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
917 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
918 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
919
920 //
22c352f8 921 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
5b33a7f5 922 Int_t bin = iPad*fMaxTime+iTimeBin;
923 allBins[iRow][bin] -= pedestalEvent;
7fef31a6 924 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
5b33a7f5 925 allBins[iRow][bin] = 0;
7fef31a6 926 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
5b33a7f5 927 allBins[iRow][bin] = 0;
22c352f8 928 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
5b33a7f5 929 allBins[iRow][bin] = 0;
930 if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
931 allBins[iRow][bin] = 0;
932 if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
22c352f8 933 }
934 }
f8aae377 935 }
22c352f8 936 }
22c352f8 937 // Now loop over rows and find clusters
938 for (fRow = 0; fRow < nRows; fRow++) {
22c352f8 939 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
ebe95b38 940 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
22c352f8 941
942 fRx = fParam->GetPadRowRadii(fSector, fRow);
943 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 944 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 945 if (fSector < kNIS)
946 fMaxPad = fParam->GetNPadsLow(fRow);
947 else
948 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 949 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 950
951 fBins = allBins[fRow];
5b33a7f5 952 fSigBins = allSigBins[fRow];
953 fNSigBins = allNSigBins[fRow];
22c352f8 954
940ed1f0 955 FindClusters(noiseROC);
ebe95b38 956 FillRow();
d101caf3 957 fRowCl->GetArray()->Clear();
22c352f8 958 nclusters += fNcluster;
959 } // End of loop to find clusters
22c352f8 960 } // End of loop over sectors
adefcaa6 961
962 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
963 delete [] allBins[iRow];
5b33a7f5 964 delete [] allSigBins[iRow];
adefcaa6 965 }
966 delete [] allBins;
5b33a7f5 967 delete [] allSigBins;
968 delete [] allNSigBins;
adefcaa6 969
9546a7ff 970 if (rawReader->GetEventId() && fOutput ){
971 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
972 }else{
973 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), nclusters);
ebe95b38 974
9546a7ff 975 }
ebe95b38 976
f8aae377 977}
978
940ed1f0 979void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 980{
940ed1f0 981
982 //
983 // add virtual charge at the edge
984 //
7041b196 985 Double_t kMaxDumpSize = 500000;
ebe95b38 986 if (!fOutput) {
987 fBDumpSignal =kFALSE;
988 }else{
989 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
990 }
991
f8aae377 992 fNcluster=0;
f8aae377 993 fLoop=1;
a1ec4d07 994 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 995 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
996 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
997 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
998 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
999 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1000 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
5b33a7f5 1001 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1002 Int_t i = fSigBins[iSig];
1003 if (i%fMaxTime<=crtime) continue;
1004 Float_t *b = &fBins[i];
940ed1f0 1005 //absolute custs
03fe1804 1006 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1007 //
1008 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
aba7fdfc 1009 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1010 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
03fe1804 1011 //
1012 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
940ed1f0 1013 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 1014 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 1015 //
1016 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
b24bd339 1017 if (noise>fRecoParam->GetMaxNoise()) continue;
940ed1f0 1018 // sigma cuts
1019 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1020 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1021 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1022
d101caf3 1023 AliTPCclusterMI c; // default cosntruction without info
f8aae377 1024 Int_t dummy=0;
1025 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 1026
f8aae377 1027 //}
1028 }
8569a2b0 1029}
65c045f0 1030
1031
940ed1f0 1032Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 1033 //
1034 // process signal on given pad - + streaming of additional information in special mode
1035 //
1036 // id[0] - sector
1037 // id[1] - row
adefcaa6 1038 // id[2] - pad
1039
1040 //
1041 // ESTIMATE pedestal and the noise
1042 //
adefcaa6 1043 const Int_t kPedMax = 100;
1044 Float_t max = 0;
1045 Float_t maxPos = 0;
1046 Int_t median = -1;
1047 Int_t count0 = 0;
1048 Int_t count1 = 0;
940ed1f0 1049 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1050 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1051 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 1052 //
adefcaa6 1053 UShort_t histo[kPedMax];
1054 memset(histo,0,kPedMax*sizeof(UShort_t));
1055 for (Int_t i=0; i<fMaxTime; i++){
1056 if (signal[i]<=0) continue;
940ed1f0 1057 if (signal[i]>max && i>firstBin) {
65c045f0 1058 max = signal[i];
1059 maxPos = i;
adefcaa6 1060 }
1061 if (signal[i]>kPedMax-1) continue;
1062 histo[int(signal[i]+0.5)]++;
1063 count0++;
65c045f0 1064 }
7fef31a6 1065 //
adefcaa6 1066 for (Int_t i=1; i<kPedMax; i++){
1067 if (count1<count0*0.5) median=i;
1068 count1+=histo[i];
1069 }
1070 // truncated mean
65c045f0 1071 //
7041b196 1072 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1073 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1074 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
adefcaa6 1075 //
1076 for (Int_t idelta=1; idelta<10; idelta++){
1077 if (median-idelta<=0) continue;
1078 if (median+idelta>kPedMax) continue;
1079 if (count06<0.6*count1){
1080 count06+=histo[median-idelta];
1081 mean06 +=histo[median-idelta]*(median-idelta);
1082 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1083 count06+=histo[median+idelta];
1084 mean06 +=histo[median+idelta]*(median+idelta);
1085 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1086 }
1087 if (count09<0.9*count1){
1088 count09+=histo[median-idelta];
1089 mean09 +=histo[median-idelta]*(median-idelta);
1090 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1091 count09+=histo[median+idelta];
1092 mean09 +=histo[median+idelta]*(median+idelta);
1093 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1094 }
adefcaa6 1095 if (count10<0.95*count1){
1096 count10+=histo[median-idelta];
1097 mean +=histo[median-idelta]*(median-idelta);
1098 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1099 count10+=histo[median+idelta];
1100 mean +=histo[median+idelta]*(median+idelta);
1101 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1102 }
1103 }
3f0af4ba 1104 if (count10) {
1105 mean /=count10;
1106 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1107 }
1108 if (count06) {
1109 mean06/=count06;
1110 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1111 }
1112 if (count09) {
1113 mean09/=count09;
1114 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1115 }
940ed1f0 1116 rmsEvent = rms09;
adefcaa6 1117 //
940ed1f0 1118 pedestalEvent = median;
adefcaa6 1119 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1120 //
1121 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1122 //
1123 // Dump mean signal info
1124 //
1125 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1126 "TimeStamp="<<fTimeStamp<<
1127 "EventType="<<fEventType<<
adefcaa6 1128 "Sector="<<uid[0]<<
1129 "Row="<<uid[1]<<
1130 "Pad="<<uid[2]<<
1131 "Max="<<max<<
1132 "MaxPos="<<maxPos<<
1133 //
1134 "Median="<<median<<
1135 "Mean="<<mean<<
1136 "RMS="<<rms<<
1137 "Mean06="<<mean06<<
1138 "RMS06="<<rms06<<
1139 "Mean09="<<mean09<<
1140 "RMS09="<<rms09<<
940ed1f0 1141 "RMSCalib="<<rmsCalib<<
1142 "PedCalib="<<pedestalCalib<<
adefcaa6 1143 "\n";
1144 //
1145 // fill pedestal histogram
1146 //
65c045f0 1147 //
adefcaa6 1148 //
1149 //
1150 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1151 Float_t *dsignal = new Float_t[nchannels];
1152 Float_t *dtime = new Float_t[nchannels];
1153 for (Int_t i=0; i<nchannels; i++){
1154 dtime[i] = i;
1155 dsignal[i] = signal[i];
1156 }
03fe1804 1157
11441748 1158 TGraph * graph=0;
7fef31a6 1159 //
a525cc34 1160 // Big signals dumping
1161 //
1162 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1163 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1164 "TimeStamp="<<fTimeStamp<<
1165 "EventType="<<fEventType<<
1166 "Sector="<<uid[0]<<
1167 "Row="<<uid[1]<<
1168 "Pad="<<uid[2]<<
1169 "Graph="<<graph<<
1170 "Max="<<max<<
1171 "MaxPos="<<maxPos<<
1172 //
1173 "Median="<<median<<
1174 "Mean="<<mean<<
1175 "RMS="<<rms<<
1176 "Mean06="<<mean06<<
1177 "RMS06="<<rms06<<
1178 "Mean09="<<mean09<<
1179 "RMS09="<<rms09<<
1180 "\n";
1181 delete graph;
7fef31a6 1182
65c045f0 1183 delete [] dsignal;
1184 delete [] dtime;
940ed1f0 1185 if (rms06>fRecoParam->GetMaxNoise()) {
1186 pedestalEvent+=1024.;
1187 return 1024+median; // sign noisy channel in debug mode
1188 }
65c045f0 1189 return median;
1190}
1191
1192
1193
03fe1804 1194
03fe1804 1195