Add cluster error and shape parameterization class + routine to fit parameters
[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;
940ed1f0 326
1c53abe2 327 return;
328 }
329 //
330 //unfolding when neccessary
331 //
332
13116aec 333 Float_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3
334 Float_t dummy[7]={0,0,0,0,0,0};
1c53abe2 335 for (Int_t di=-3;di<=3;di++){
336 matrix2[di+3] = &bins[k+di*max];
337 if ((k+di*max)<3) matrix2[di+3] = &dummy[3];
338 if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3];
339 }
340 Float_t vmatrix2[5][5];
341 Float_t sumu;
342 Float_t overlap;
343 UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
344 //
345 // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
346 meani +=i0;
347 meanj +=j0;
348 //set cluster parameters
349 c.SetQ(sumu);
350 c.SetY(meani*fPadWidth);
351 c.SetZ(meanj*fZWidth);
352 c.SetSigmaY2(mi2);
353 c.SetSigmaZ2(mj2);
354 c.SetType(Char_t(overlap)+1);
355 AddCluster(c);
356
357 //unfolding 2
358 meani-=i0;
359 meanj-=j0;
360 if (gDebug>4)
361 printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
362}
363
364
365
13116aec 366void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
1c53abe2 367 Float_t & sumu, Float_t & overlap )
368{
369 //
370 //unfold cluster from input matrix
371 //data corresponding to cluster writen in recmatrix
372 //output meani and meanj
373
374 //take separatelly y and z
375
376 Float_t sum3i[7] = {0,0,0,0,0,0,0};
377 Float_t sum3j[7] = {0,0,0,0,0,0,0};
378
379 for (Int_t k =0;k<7;k++)
380 for (Int_t l = -1; l<=1;l++){
381 sum3i[k]+=matrix2[k][l];
382 sum3j[k]+=matrix2[l+3][k-3];
383 }
384 Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
385 //
386 //unfold y
387 Float_t sum3wi = 0; //charge minus overlap
388 Float_t sum3wio = 0; //full charge
389 Float_t sum3iw = 0; //sum for mean value
390 for (Int_t dk=-1;dk<=1;dk++){
391 sum3wio+=sum3i[dk+3];
392 if (dk==0){
393 sum3wi+=sum3i[dk+3];
394 }
395 else{
396 Float_t ratio =1;
397 if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
398 sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
399 Float_t xm2 = sum3i[-dk+3];
400 Float_t xm1 = sum3i[+3];
401 Float_t x1 = sum3i[2*dk+3];
402 Float_t x2 = sum3i[3*dk+3];
cc5e9db0 403 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
404 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 405 ratio = w11/(w11+w12);
406 for (Int_t dl=-1;dl<=1;dl++)
407 mratio[dk+1][dl+1] *= ratio;
408 }
409 Float_t amp = sum3i[dk+3]*ratio;
410 sum3wi+=amp;
411 sum3iw+= dk*amp;
412 }
413 }
414 meani = sum3iw/sum3wi;
415 Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
416
417
418
419 //unfold z
420 Float_t sum3wj = 0; //charge minus overlap
421 Float_t sum3wjo = 0; //full charge
422 Float_t sum3jw = 0; //sum for mean value
423 for (Int_t dk=-1;dk<=1;dk++){
424 sum3wjo+=sum3j[dk+3];
425 if (dk==0){
426 sum3wj+=sum3j[dk+3];
427 }
428 else{
429 Float_t ratio =1;
430 if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
431 (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
432 Float_t xm2 = sum3j[-dk+3];
433 Float_t xm1 = sum3j[+3];
434 Float_t x1 = sum3j[2*dk+3];
435 Float_t x2 = sum3j[3*dk+3];
cc5e9db0 436 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
437 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
1c53abe2 438 ratio = w11/(w11+w12);
439 for (Int_t dl=-1;dl<=1;dl++)
440 mratio[dl+1][dk+1] *= ratio;
441 }
442 Float_t amp = sum3j[dk+3]*ratio;
443 sum3wj+=amp;
444 sum3jw+= dk*amp;
445 }
446 }
447 meanj = sum3jw/sum3wj;
448 Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
449 overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
450 sumu = (sum3wj+sum3wi)/2.;
451
452 if (overlap ==3) {
453 //if not overlap detected remove everything
454 for (Int_t di =-2; di<=2;di++)
455 for (Int_t dj =-2; dj<=2;dj++){
456 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
457 }
458 }
459 else{
460 for (Int_t di =-1; di<=1;di++)
461 for (Int_t dj =-1; dj<=1;dj++){
462 Float_t ratio =1;
463 if (mratio[di+1][dj+1]==1){
464 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
465 if (TMath::Abs(di)+TMath::Abs(dj)>1){
466 recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
467 recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
468 }
469 recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
470 }
471 else
472 {
473 //if we have overlap in direction
474 recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
475 if (TMath::Abs(di)+TMath::Abs(dj)>1){
cc5e9db0 476 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
1c53abe2 477 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
478 //
cc5e9db0 479 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
1c53abe2 480 recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
481 }
482 else{
483 ratio = recmatrix[di+2][dj+2]/matrix2[3][0];
484 recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
485 }
486 }
487 }
488 }
489 if (gDebug>4)
490 printf("%f\n", recmatrix[2][2]);
491
492}
493
494Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
495{
496 //
497 // estimate max
498 Float_t sumteor= 0;
499 Float_t sumamp = 0;
500
501 for (Int_t di = -1;di<=1;di++)
502 for (Int_t dj = -1;dj<=1;dj++){
503 if (vmatrix[2+di][2+dj]>2){
504 Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
505 sumteor += teor*vmatrix[2+di][2+dj];
506 sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
507 }
508 }
509 Float_t max = sumamp/sumteor;
510 return max;
511}
512
513void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
514 //
515 // transform cluster to the global coordinata
516 // add the cluster to the array
517 //
518 Float_t meani = c.GetY()/fPadWidth;
519 Float_t meanj = c.GetZ()/fZWidth;
520
521 Int_t ki = TMath::Nint(meani-3);
522 if (ki<0) ki=0;
523 if (ki>=fMaxPad) ki = fMaxPad-1;
524 Int_t kj = TMath::Nint(meanj-3);
525 if (kj<0) kj=0;
526 if (kj>=fMaxTime-3) kj=fMaxTime-4;
527 // ki and kj shifted to "real" coordinata
f8aae377 528 if (fRowDig) {
529 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
530 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
531 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
532 }
1c53abe2 533
534
535 Float_t s2 = c.GetSigmaY2();
536 Float_t w=fParam->GetPadPitchWidth(fSector);
537
538 c.SetSigmaY2(s2*w*w);
539 s2 = c.GetSigmaZ2();
540 w=fZWidth;
541 c.SetSigmaZ2(s2*w*w);
542 c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
adefcaa6 543 if (!fRecoParam->GetBYMirror()){
544 if (fSector%36>17){
545 c.SetY(-(meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
546 }
547 }
1c53abe2 548 c.SetZ(fZWidth*(meanj-3));
753797ce 549 c.SetZ(c.GetZ() - 3.*fParam->GetZSigma() + fParam->GetNTBinsL1()*fParam->GetZWidth()); // PASA delay + L1 delay
1c53abe2 550 c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
508541c7 551 c.SetX(fRx);
552 c.SetDetector(fSector);
553 c.SetRow(fRow);
554
1c53abe2 555 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
556 //c.SetSigmaY2(c.GetSigmaY2()*25.);
557 //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
558 c.SetType(-(c.GetType()+3)); //edge clusters
559 }
560 if (fLoop==2) c.SetType(100);
561
562 TClonesArray * arr = fRowCl->GetArray();
45bcf167 563 // AliTPCclusterMI * cl =
564 new ((*arr)[fNcluster]) AliTPCclusterMI(c);
1c53abe2 565
566 fNcluster++;
567}
568
569
570//_____________________________________________________________________________
f8aae377 571void AliTPCclustererMI::Digits2Clusters()
1c53abe2 572{
573 //-----------------------------------------------------------------
574 // This is a simple cluster finder.
575 //-----------------------------------------------------------------
1c53abe2 576
f8aae377 577 if (!fInput) {
578 Error("Digits2Clusters", "input tree not initialised");
1c53abe2 579 return;
580 }
581
f8aae377 582 if (!fOutput) {
583 Error("Digits2Clusters", "output tree not initialised");
584 return;
1c53abe2 585 }
586
13116aec 587 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 588 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
13116aec 589
1c53abe2 590 AliSimDigits digarr, *dummy=&digarr;
591 fRowDig = dummy;
592 fInput->GetBranch("Segment")->SetAddress(&dummy);
593 Stat_t nentries = fInput->GetEntries();
594
f8aae377 595 fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after
1c53abe2 596
1c53abe2 597 Int_t nclusters = 0;
13116aec 598
1c53abe2 599 for (Int_t n=0; n<nentries; n++) {
600 fInput->GetEvent(n);
508541c7 601 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
1c53abe2 602 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
603 continue;
604 }
508541c7 605 Int_t row = fRow;
13116aec 606 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
940ed1f0 607 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
13116aec 608 //
1c53abe2 609 AliTPCClustersRow *clrow= new AliTPCClustersRow();
610 fRowCl = clrow;
611 clrow->SetClass("AliTPCclusterMI");
612 clrow->SetArray(1);
613
614 clrow->SetID(digarr.GetID());
615 fOutput->GetBranch("Segment")->SetAddress(&clrow);
f8aae377 616 fRx=fParam->GetPadRowRadii(fSector,row);
1c53abe2 617
618
f8aae377 619 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
1c53abe2 620 fZWidth = fParam->GetZWidth();
621 if (fSector < kNIS) {
f8aae377 622 fMaxPad = fParam->GetNPadsLow(row);
1c53abe2 623 fSign = (fSector < kNIS/2) ? 1 : -1;
f8aae377 624 fPadLength = fParam->GetPadPitchLength(fSector,row);
625 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 626 } else {
f8aae377 627 fMaxPad = fParam->GetNPadsUp(row);
1c53abe2 628 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
f8aae377 629 fPadLength = fParam->GetPadPitchLength(fSector,row);
630 fPadWidth = fParam->GetPadPitchWidth();
1c53abe2 631 }
632
633
634 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
13116aec 635 fBins =new Float_t[fMaxBin];
636 fResBins =new Float_t[fMaxBin]; //fBins with residuals after 1 finder loop
637 memset(fBins,0,sizeof(Float_t)*fMaxBin);
638 memset(fResBins,0,sizeof(Float_t)*fMaxBin);
1c53abe2 639
640 if (digarr.First()) //MI change
641 do {
13116aec 642 Float_t dig=digarr.CurrentDigit();
f8aae377 643 if (dig<=fParam->GetZeroSup()) continue;
1c53abe2 644 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
9d4f75a9 645 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
13116aec 646 fBins[i*fMaxTime+j]=dig/gain;
1c53abe2 647 } while (digarr.Next());
648 digarr.ExpandTrackBuffer();
649
940ed1f0 650 FindClusters(noiseROC);
8569a2b0 651
652 fOutput->Fill();
88cb7938 653 delete clrow;
654 nclusters+=fNcluster;
8569a2b0 655 delete[] fBins;
656 delete[] fResBins;
88cb7938 657 }
f8aae377 658
19dd5b2f 659 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
f8aae377 660}
661
662void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
663{
664//-----------------------------------------------------------------
22c352f8 665// This is a cluster finder for the TPC raw data.
666// The method assumes NO ordering of the altro channels.
667// The pedestal subtraction can be switched on and off
668// using an option of the TPC reconstructor
f8aae377 669//-----------------------------------------------------------------
670
671 if (!fOutput) {
672 Error("Digits2Clusters", "output tree not initialised");
673 return;
674 }
675
f8aae377 676 fRowDig = NULL;
adefcaa6 677 AliTPCROC * roc = AliTPCROC::Instance();
22c352f8 678 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
940ed1f0 679 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
680 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
adefcaa6 681 AliTPCRawStream input(rawReader);
ef4ad662 682 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
683 if (fEventHeader){
684 fTimeStamp = fEventHeader->Get("Timestamp");
685 fEventType = fEventHeader->Get("Type");
686 }
687
22c352f8 688
f8aae377 689 Int_t nclusters = 0;
88cb7938 690
f8aae377 691 fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after
692 const Int_t kNIS = fParam->GetNInnerSector();
693 const Int_t kNOS = fParam->GetNOuterSector();
694 const Int_t kNS = kNIS + kNOS;
695 fZWidth = fParam->GetZWidth();
696 Int_t zeroSup = fParam->GetZeroSup();
adefcaa6 697 //
698 //alocate memory for sector - maximal case
699 //
22c352f8 700 Float_t** allBins = NULL;
701 Float_t** allBinsRes = NULL;
adefcaa6 702 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
703 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
704 allBins = new Float_t*[nRowsMax];
705 allBinsRes = new Float_t*[nRowsMax];
706 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
707 //
708 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
709 allBins[iRow] = new Float_t[maxBin];
710 allBinsRes[iRow] = new Float_t[maxBin];
711 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
adefcaa6 712 }
713 //
22c352f8 714 // Loop over sectors
adefcaa6 715 //
22c352f8 716 for(fSector = 0; fSector < kNS; fSector++) {
f8aae377 717
940ed1f0 718 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
719 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
720 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
22c352f8 721
722 Int_t nRows = 0;
723 Int_t nDDLs = 0, indexDDL = 0;
724 if (fSector < kNIS) {
725 nRows = fParam->GetNRowLow();
726 fSign = (fSector < kNIS/2) ? 1 : -1;
727 nDDLs = 2;
728 indexDDL = fSector * 2;
729 }
730 else {
731 nRows = fParam->GetNRowUp();
732 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
733 nDDLs = 4;
734 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
735 }
736
22c352f8 737 for (Int_t iRow = 0; iRow < nRows; iRow++) {
738 Int_t maxPad;
739 if (fSector < kNIS)
740 maxPad = fParam->GetNPadsLow(iRow);
741 else
742 maxPad = fParam->GetNPadsUp(iRow);
743
744 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
22c352f8 745 memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
22c352f8 746 }
747
748 // Loas the raw data for corresponding DDLs
749 rawReader->Reset();
22c352f8 750 input.SetOldRCUFormat(fIsOldRCUFormat);
362c9d61 751 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
adefcaa6 752 Int_t digCounter=0;
22c352f8 753 // Begin loop over altro data
adefcaa6 754 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
755 Float_t gain =1;
756 Int_t lastPad=-1;
22c352f8 757 while (input.Next()) {
adefcaa6 758 digCounter++;
22c352f8 759 if (input.GetSector() != fSector)
760 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
761
22c352f8 762
763 Int_t iRow = input.GetRow();
764 if (iRow < 0 || iRow >= nRows)
765 AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
766 iRow, 0, nRows -1));
adefcaa6 767 //pad
768 Int_t iPad = input.GetPad();
769 if (iPad < 0 || iPad >= nPadsMax)
22c352f8 770 AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
adefcaa6 771 iPad, 0, nPadsMax-1));
772 if (iPad!=lastPad){
773 gain = gainROC->GetValue(iRow,iPad);
774 lastPad = iPad;
775 }
776 iPad+=3;
777 //time
778 Int_t iTimeBin = input.GetTime();
779 if ( iTimeBin < 0 || iTimeBin >= fParam->GetMaxTBin())
22c352f8 780 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
adefcaa6 781 iTimeBin, 0, iTimeBin -1));
782 iTimeBin+=3;
783 //signal
22c352f8 784 Float_t signal = input.GetSignal();
adefcaa6 785 if (!calcPedestal && signal <= zeroSup) continue;
940ed1f0 786 if (!calcPedestal) {
787 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal/gain;
788 }else{
789 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
790 }
791 allBins[iRow][iPad*fMaxTime+0]=1.; // pad with signal
22c352f8 792 } // End of the loop over altro data
adefcaa6 793 //
794 //
22c352f8 795 // Now loop over rows and perform pedestal subtraction
adefcaa6 796 if (digCounter==0) continue;
194b0609 797 // if (fPedSubtraction) {
adefcaa6 798 if (calcPedestal) {
22c352f8 799 for (Int_t iRow = 0; iRow < nRows; iRow++) {
800 Int_t maxPad;
801 if (fSector < kNIS)
802 maxPad = fParam->GetNPadsLow(iRow);
803 else
804 maxPad = fParam->GetNPadsUp(iRow);
805
940ed1f0 806 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
807 if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
22c352f8 808 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
65c045f0 809 //Float_t pedestal = TMath::Median(fMaxTime, p);
810 Int_t id[3] = {fSector, iRow, iPad-3};
940ed1f0 811 // calib values
812 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
813 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
814 Double_t rmsEvent = rmsCalib;
815 Double_t pedestalEvent = pedestalCalib;
816 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
817 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
818 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
819
820 //
22c352f8 821 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
940ed1f0 822 allBins[iRow][iPad*fMaxTime+iTimeBin] -= pedestalEvent;
7fef31a6 823 if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
824 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
825 if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
826 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
22c352f8 827 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
828 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
940ed1f0 829 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < 3.0*rmsEvent) // 3 sigma cut on RMS
ef4ad662 830 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
22c352f8 831 }
832 }
f8aae377 833 }
22c352f8 834 }
22c352f8 835 // Now loop over rows and find clusters
836 for (fRow = 0; fRow < nRows; fRow++) {
837 fRowCl = new AliTPCClustersRow;
838 fRowCl->SetClass("AliTPCclusterMI");
839 fRowCl->SetArray(1);
840 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
841 fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
842
843 fRx = fParam->GetPadRowRadii(fSector, fRow);
844 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
f8aae377 845 fPadWidth = fParam->GetPadPitchWidth();
22c352f8 846 if (fSector < kNIS)
847 fMaxPad = fParam->GetNPadsLow(fRow);
848 else
849 fMaxPad = fParam->GetNPadsUp(fRow);
f8aae377 850 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
22c352f8 851
852 fBins = allBins[fRow];
853 fResBins = allBinsRes[fRow];
854
940ed1f0 855 FindClusters(noiseROC);
22c352f8 856
857 fOutput->Fill();
858 delete fRowCl;
859 nclusters += fNcluster;
860 } // End of loop to find clusters
22c352f8 861 } // End of loop over sectors
adefcaa6 862
863 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
864 delete [] allBins[iRow];
865 delete [] allBinsRes[iRow];
866 }
867 delete [] allBins;
868 delete [] allBinsRes;
869
940ed1f0 870 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
22c352f8 871
f8aae377 872}
873
940ed1f0 874void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
f8aae377 875{
940ed1f0 876
877 //
878 // add virtual charge at the edge
879 //
880 if (0) for (Int_t i=0; i<fMaxTime; i++){
f8aae377 881 Float_t amp1 = fBins[i+3*fMaxTime];
882 Float_t amp0 =0;
883 if (amp1>0){
884 Float_t amp2 = fBins[i+4*fMaxTime];
885 if (amp2==0) amp2=0.5;
886 Float_t sigma2 = GetSigmaY2(i);
887 amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
888 if (gDebug>4) printf("\n%f\n",amp0);
889 }
13116aec 890 fBins[i+2*fMaxTime] = amp0;
f8aae377 891 amp0 = 0;
892 amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
893 if (amp1>0){
894 Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
895 if (amp2==0) amp2=0.5;
896 Float_t sigma2 = GetSigmaY2(i);
897 amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
898 if (gDebug>4) printf("\n%f\n",amp0);
899 }
13116aec 900 fBins[(fMaxPad+3)*fMaxTime+i] = amp0;
f8aae377 901 }
940ed1f0 902 memcpy(fResBins,fBins, fMaxBin*sizeof(Float_t));
903 //
904 //
f8aae377 905 //
906 fNcluster=0;
f8aae377 907 fLoop=1;
13116aec 908 Float_t *b=&fBins[-1]+2*fMaxTime;
194b0609 909 Int_t crtime = Int_t((fParam->GetZLength()-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
940ed1f0 910 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
911 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
912 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
913 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
914 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
915 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
f8aae377 916 for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
917 b++;
f8aae377 918 if (i%fMaxTime<crtime) {
919 Int_t delta = -(i%fMaxTime)+crtime;
920 b+=delta;
921 i+=delta;
922 continue;
923 }
940ed1f0 924 //absolute custs
925 if (b[0]<minMaxCutAbs) continue; //threshold form maxima
926 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
927 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
928 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up town TRF
929 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
f8aae377 930 if (!IsMaximum(*b,fMaxTime,b)) continue;
940ed1f0 931 //
932 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
933 // sigma cuts
934 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
935 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
936 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
937
f8aae377 938 AliTPCclusterMI c;
939 Int_t dummy=0;
940 MakeCluster(i, fMaxTime, fBins, dummy,c);
940ed1f0 941
f8aae377 942 //}
943 }
8569a2b0 944}
65c045f0 945
946
940ed1f0 947Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
65c045f0 948 //
949 // process signal on given pad - + streaming of additional information in special mode
950 //
951 // id[0] - sector
952 // id[1] - row
adefcaa6 953 // id[2] - pad
954
955 //
956 // ESTIMATE pedestal and the noise
957 //
65c045f0 958 Int_t offset =100;
adefcaa6 959 const Int_t kPedMax = 100;
960 Float_t max = 0;
961 Float_t maxPos = 0;
962 Int_t median = -1;
963 Int_t count0 = 0;
964 Int_t count1 = 0;
940ed1f0 965 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
966 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
967 Int_t firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
65c045f0 968 //
adefcaa6 969 UShort_t histo[kPedMax];
970 memset(histo,0,kPedMax*sizeof(UShort_t));
971 for (Int_t i=0; i<fMaxTime; i++){
972 if (signal[i]<=0) continue;
940ed1f0 973 if (signal[i]>max && i>firstBin) {
65c045f0 974 max = signal[i];
975 maxPos = i;
adefcaa6 976 }
977 if (signal[i]>kPedMax-1) continue;
978 histo[int(signal[i]+0.5)]++;
979 count0++;
65c045f0 980 }
7fef31a6 981 //
adefcaa6 982 for (Int_t i=1; i<kPedMax; i++){
983 if (count1<count0*0.5) median=i;
984 count1+=histo[i];
985 }
986 // truncated mean
65c045f0 987 //
adefcaa6 988 Double_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
989 Double_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
990 Double_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
991 //
992 for (Int_t idelta=1; idelta<10; idelta++){
993 if (median-idelta<=0) continue;
994 if (median+idelta>kPedMax) continue;
995 if (count06<0.6*count1){
996 count06+=histo[median-idelta];
997 mean06 +=histo[median-idelta]*(median-idelta);
998 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
999 count06+=histo[median+idelta];
1000 mean06 +=histo[median+idelta]*(median+idelta);
1001 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1002 }
1003 if (count09<0.9*count1){
1004 count09+=histo[median-idelta];
1005 mean09 +=histo[median-idelta]*(median-idelta);
1006 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1007 count09+=histo[median+idelta];
1008 mean09 +=histo[median+idelta]*(median+idelta);
1009 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
65c045f0 1010 }
adefcaa6 1011 if (count10<0.95*count1){
1012 count10+=histo[median-idelta];
1013 mean +=histo[median-idelta]*(median-idelta);
1014 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1015 count10+=histo[median+idelta];
1016 mean +=histo[median+idelta]*(median+idelta);
1017 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1018 }
1019 }
1020 mean /=count10;
1021 mean06/=count06;
1022 mean09/=count09;
1023 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1024 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
940ed1f0 1025 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1026 rmsEvent = rms09;
adefcaa6 1027 //
940ed1f0 1028 pedestalEvent = median;
adefcaa6 1029 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1030 //
1031 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1032 //
1033 // Dump mean signal info
1034 //
1035 (*fDebugStreamer)<<"Signal"<<
ef4ad662 1036 "TimeStamp="<<fTimeStamp<<
1037 "EventType="<<fEventType<<
adefcaa6 1038 "Sector="<<uid[0]<<
1039 "Row="<<uid[1]<<
1040 "Pad="<<uid[2]<<
1041 "Max="<<max<<
1042 "MaxPos="<<maxPos<<
1043 //
1044 "Median="<<median<<
1045 "Mean="<<mean<<
1046 "RMS="<<rms<<
1047 "Mean06="<<mean06<<
1048 "RMS06="<<rms06<<
1049 "Mean09="<<mean09<<
1050 "RMS09="<<rms09<<
940ed1f0 1051 "RMSCalib="<<rmsCalib<<
1052 "PedCalib="<<pedestalCalib<<
adefcaa6 1053 "\n";
1054 //
1055 // fill pedestal histogram
1056 //
1057 AliTPCROC * roc = AliTPCROC::Instance();
1058 if (!fAmplitudeHisto){
1059 fAmplitudeHisto = new TObjArray(72);
1060 }
1061 //
1062 if (uid[0]<roc->GetNSectors()
1063 && uid[1]< roc->GetNRows(uid[0]) &&
1064 uid[2] <roc->GetNPads(uid[0], uid[1])){
65c045f0 1065 TObjArray * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
1066 if (!sectorArray){
adefcaa6 1067 Int_t npads =roc->GetNChannels(uid[0]);
65c045f0 1068 sectorArray = new TObjArray(npads);
1069 fAmplitudeHisto->AddAt(sectorArray, uid[0]);
1070 }
adefcaa6 1071 Int_t position = uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
65c045f0 1072 TH1F * histo = (TH1F*)sectorArray->UncheckedAt(position);
1073 if (!histo){
1074 char hname[100];
1075 sprintf(hname,"Amp_%d_%d_%d",uid[0],uid[1],uid[2]);
1076 TFile * backup = gFile;
1077 fDebugStreamer->GetFile()->cd();
194b0609 1078 histo = new TH1F(hname, hname, 100, 5,100);
65c045f0 1079 //histo->SetDirectory(0); // histogram not connected to directory -(File)
1080 sectorArray->AddAt(histo, position);
1081 if (backup) backup->cd();
1082 }
1083 for (Int_t i=0; i<nchannels; i++){
adefcaa6 1084 histo->Fill(signal[i]);
65c045f0 1085 }
1086 }
1087 //
adefcaa6 1088 //
1089 //
1090 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1091 Float_t *dsignal = new Float_t[nchannels];
1092 Float_t *dtime = new Float_t[nchannels];
1093 for (Int_t i=0; i<nchannels; i++){
1094 dtime[i] = i;
1095 dsignal[i] = signal[i];
1096 }
1097 //
1098 //
65c045f0 1099 TGraph * graph;
1100 Bool_t random = (gRandom->Rndm()<0.0001);
1101 if (max-median>kMin || rms06>2.*fParam->GetZeroSup() || random){
1102 graph =new TGraph(nchannels, dtime, dsignal);
1103 if (rms06>2.*fParam->GetZeroSup() || random)
1104 (*fDebugStreamer)<<"SignalN"<< //noise pads - or random sample of pads
ef4ad662 1105 "TimeStamp="<<fTimeStamp<<
1106 "EventType="<<fEventType<<
65c045f0 1107 "Sector="<<uid[0]<<
1108 "Row="<<uid[1]<<
1109 "Pad="<<uid[2]<<
864c1803 1110 "Graph="<<graph<<
65c045f0 1111 "Max="<<max<<
1112 "MaxPos="<<maxPos<<
1113 //
1114 "Median="<<median<<
1115 "Mean="<<mean<<
1116 "RMS="<<rms<<
1117 "Mean06="<<mean06<<
1118 "RMS06="<<rms06<<
65c045f0 1119 "Mean09="<<mean09<<
1120 "RMS09="<<rms09<<
1121 "\n";
fe5055e5 1122 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
65c045f0 1123 (*fDebugStreamer)<<"SignalB"<< // pads with signal
ef4ad662 1124 "TimeStamp="<<fTimeStamp<<
1125 "EventType="<<fEventType<<
65c045f0 1126 "Sector="<<uid[0]<<
1127 "Row="<<uid[1]<<
1128 "Pad="<<uid[2]<<
864c1803 1129 "Graph="<<graph<<
65c045f0 1130 "Max="<<max<<
1131 "MaxPos="<<maxPos<<
1132 //
1133 "Median="<<median<<
1134 "Mean="<<mean<<
1135 "RMS="<<rms<<
1136 "Mean06="<<mean06<<
1137 "RMS06="<<rms06<<
65c045f0 1138 "Mean09="<<mean09<<
1139 "RMS09="<<rms09<<
1140 "\n";
1141 delete graph;
1142 }
1143
7fef31a6 1144 //
1145 //
1146 // Central Electrode signal analysis
1147 //
1148 Double_t ceQmax =0, ceQsum=0, ceTime=0;
1149 Double_t cemean = mean06, cerms=rms06 ;
1150 Int_t cemaxpos= 0;
1151 Double_t ceThreshold=5.*cerms;
1152 Double_t ceSumThreshold=8.*cerms;
97f127bb 1153 const Int_t kCemin=5; // range for the analysis of the ce signal +- channels from the peak
1154 const Int_t kCemax=5;
adefcaa6 1155 for (Int_t i=nchannels-2; i>nchannels/2; i--){
7fef31a6 1156 if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){
1157 cemaxpos=i;
1158 break;
1159 }
1160 }
1161 if (cemaxpos!=0){
adefcaa6 1162 ceQmax = 0;
1163 Int_t cemaxpos2=0;
1164 for (Int_t i=cemaxpos-20; i<cemaxpos+5; i++){
1165 if (i<0 || i>nchannels-1) continue;
1166 Double_t val=dsignal[i]- cemean;
1167 if (val>ceQmax){
1168 cemaxpos2=i;
1169 ceQmax = val;
7fef31a6 1170 }
adefcaa6 1171 }
1172 cemaxpos = cemaxpos2;
1173
1174 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
1175 if (i>0 && i<nchannels&&dsignal[i]- cemean>0){
1176 Double_t val=dsignal[i]- cemean;
1177 ceTime+=val*dtime[i];
1178 ceQsum+=val;
1179 if (val>ceQmax) ceQmax=val;
7fef31a6 1180 }
adefcaa6 1181 }
1182 if (ceQmax&&ceQsum>ceSumThreshold) {
1183 ceTime/=ceQsum;
1184 (*fDebugStreamer)<<"Signalce"<<
ef4ad662 1185 "TimeStamp="<<fTimeStamp<<
1186 "EventType="<<fEventType<<
adefcaa6 1187 "Sector="<<uid[0]<<
1188 "Row="<<uid[1]<<
1189 "Pad="<<uid[2]<<
1190 "Max="<<ceQmax<<
1191 "Qsum="<<ceQsum<<
1192 "Time="<<ceTime<<
1193 "RMS06="<<rms06<<
1194 //
1195 "\n";
1196 }
7fef31a6 1197 }
1198 // end of ce signal analysis
1199 //
1200
1201 //
1202 // Gating grid signal analysis
1203 //
1204 Double_t ggQmax =0, ggQsum=0, ggTime=0;
1205 Double_t ggmean = mean06, ggrms=rms06 ;
1206 Int_t ggmaxpos= 0;
1207 Double_t ggThreshold=5.*ggrms;
1208 Double_t ggSumThreshold=8.*ggrms;
bd85a5e8 1209
adefcaa6 1210 for (Int_t i=1; i<nchannels/4; i++){
bd85a5e8 1211 if ( (dsignal[i]-mean06)>ggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] &&
1212 (dsignal[i]+dsignal[i+1]+dsignal[i-1]-3*mean06)>ggSumThreshold){
7fef31a6 1213 ggmaxpos=i;
bd85a5e8 1214 if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
7fef31a6 1215 break;
1216 }
1217 }
1218 if (ggmaxpos!=0){
bd85a5e8 1219 for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){
1220 if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
7fef31a6 1221 Double_t val=dsignal[i]- ggmean;
1222 ggTime+=val*dtime[i];
1223 ggQsum+=val;
1224 if (val>ggQmax) ggQmax=val;
1225 }
1226 }
1227 if (ggQmax&&ggQsum>ggSumThreshold) {
1228 ggTime/=ggQsum;
1229 (*fDebugStreamer)<<"Signalgg"<<
ef4ad662 1230 "TimeStamp="<<fTimeStamp<<
1231 "EventType="<<fEventType<<
7fef31a6 1232 "Sector="<<uid[0]<<
1233 "Row="<<uid[1]<<
1234 "Pad="<<uid[2]<<
1235 "Max="<<ggQmax<<
1236 "Qsum="<<ggQsum<<
1237 "Time="<<ggTime<<
bd85a5e8 1238 "RMS06="<<rms06<<
7fef31a6 1239 //
1240 "\n";
1241 }
1242 }
1243 // end of gg signal analysis
1244
1245
65c045f0 1246 delete [] dsignal;
1247 delete [] dtime;
940ed1f0 1248 if (rms06>fRecoParam->GetMaxNoise()) {
1249 pedestalEvent+=1024.;
1250 return 1024+median; // sign noisy channel in debug mode
1251 }
65c045f0 1252 return median;
1253}
1254
1255
1256
1257void AliTPCclustererMI::DumpHistos(){
adefcaa6 1258 //
1259 // Dump histogram information
1260 //
65c045f0 1261 if (!fAmplitudeHisto) return;
adefcaa6 1262 AliTPCROC * roc = AliTPCROC::Instance();
65c045f0 1263 for (UInt_t isector=0; isector<AliTPCROC::Instance()->GetNSectors(); isector++){
1264 TObjArray * array = (TObjArray*)fAmplitudeHisto->UncheckedAt(isector);
1265 if (!array) continue;
1266 for (UInt_t ipad = 0; ipad <(UInt_t)array->GetEntriesFast(); ipad++){
1267 TH1F * histo = (TH1F*) array->UncheckedAt(ipad);
1268 if (!histo) continue;
1269 if (histo->GetEntries()<100) continue;
1270 histo->Fit("gaus","q");
1271 Float_t mean = histo->GetMean();
1272 Float_t rms = histo->GetRMS();
1273 Float_t gmean = histo->GetFunction("gaus")->GetParameter(1);
1274 Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2);
1275 Float_t max = histo->GetFunction("gaus")->GetParameter(0);
1276
1277 // get pad number
1278 UInt_t row=0, pad =0;
adefcaa6 1279 const UInt_t *indexes =roc->GetRowIndexes(isector);
1280 for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
940ed1f0 1281 if (indexes[irow]<=ipad){
65c045f0 1282 row = irow;
1283 pad = ipad-indexes[irow];
1284 }
1285 }
1286 Int_t rpad = pad - (AliTPCROC::Instance()->GetNPads(isector,row))/2;
1287 //
1288 (*fDebugStreamer)<<"Fit"<<
ef4ad662 1289 "TimeStamp="<<fTimeStamp<<
1290 "EventType="<<fEventType<<
65c045f0 1291 "Sector="<<isector<<
1292 "Row="<<row<<
1293 "Pad="<<pad<<
1294 "RPad="<<rpad<<
1295 "Max="<<max<<
1296 "Mean="<<mean<<
1297 "RMS="<<rms<<
1298 "GMean="<<gmean<<
1299 "GSigma="<<gsigma<<
1300 "\n";
1301 if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
1302 }
1303 }
1304}