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