]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TPC/AliTPCclustererMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
AliHMPIDDigitN no longer needed
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//-------------------------------------------------------
19// Implementation of the TPC clusterer
20//
21// Origin: Marian Ivanov
22//-------------------------------------------------------
23
24#include "AliTPCReconstructor.h"
25#include "AliTPCclustererMI.h"
26#include "AliTPCclusterMI.h"
27#include <TObjArray.h>
28#include <TFile.h>
29#include "TGraph.h"
30#include "TF1.h"
31#include "TRandom.h"
32#include "AliMathBase.h"
33
34#include "AliTPCClustersArray.h"
35#include "AliTPCClustersRow.h"
36#include "AliDigits.h"
37#include "AliSimDigits.h"
38#include "AliTPCParam.h"
39#include "AliTPCRecoParam.h"
40#include "AliRawReader.h"
41#include "AliTPCRawStream.h"
42#include "AliRawEventHeaderBase.h"
43#include "AliRunLoader.h"
44#include "AliLoader.h"
45#include "Riostream.h"
46#include <TTree.h>
47#include "AliTPCcalibDB.h"
48#include "AliTPCCalPad.h"
49#include "AliTPCCalROC.h"
50#include "TTreeStream.h"
51#include "AliLog.h"
52
53
54ClassImp(AliTPCclustererMI)
55
56
57
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),
74 fEventHeader(0),
75 fTimeStamp(0),
76 fEventType(0),
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)
86{
87 //
88 // COSNTRUCTOR
89 // param - tpc parameters for given file
90 // recoparam - reconstruction parameters
91 //
92 fIsOldRCUFormat = kFALSE;
93 fInput =0;
94 fOutput=0;
95 fParam = par;
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 }
103 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
104 fAmplitudeHisto = 0;
105}
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//______________________________________________________________
152AliTPCclustererMI::~AliTPCclustererMI(){
153 DumpHistos();
154 if (fAmplitudeHisto) delete fAmplitudeHisto;
155 if (fDebugStreamer) delete fDebugStreamer;
156}
157
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
186 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
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
197 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
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
208void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
209AliTPCclusterMI &c)
210{
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 //
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};
222 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
223 Float_t * resmatrix[5];
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);
235 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
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;
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
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++){
322 resmatrix[di+2][dj] -= vmatrix[di+2][dj+2];
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
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};
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
365void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
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];
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.);
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];
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.);
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){
475 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
476 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
477 //
478 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
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
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 }
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));
542 if (!fRecoParam->GetBYMirror()){
543 if (fSector%36>17){
544 c.SetY(-(meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
545 }
546 }
547 c.SetZ(fZWidth*(meanj-3));
548 c.SetZ(c.GetZ() - 3.*fParam->GetZSigma() + fParam->GetNTBinsL1()*fParam->GetZWidth()); // PASA delay + L1 delay
549 c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
550 c.SetX(fRx);
551 c.SetDetector(fSector);
552 c.SetRow(fRow);
553
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();
562 // AliTPCclusterMI * cl =
563 new ((*arr)[fNcluster]) AliTPCclusterMI(c);
564
565 fNcluster++;
566}
567
568
569//_____________________________________________________________________________
570void AliTPCclustererMI::Digits2Clusters()
571{
572 //-----------------------------------------------------------------
573 // This is a simple cluster finder.
574 //-----------------------------------------------------------------
575
576 if (!fInput) {
577 Error("Digits2Clusters", "input tree not initialised");
578 return;
579 }
580
581 if (!fOutput) {
582 Error("Digits2Clusters", "output tree not initialised");
583 return;
584 }
585
586 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
587
588 AliSimDigits digarr, *dummy=&digarr;
589 fRowDig = dummy;
590 fInput->GetBranch("Segment")->SetAddress(&dummy);
591 Stat_t nentries = fInput->GetEntries();
592
593 fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after
594
595 Int_t nclusters = 0;
596
597 for (Int_t n=0; n<nentries; n++) {
598 fInput->GetEvent(n);
599 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
600 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
601 continue;
602 }
603 Int_t row = fRow;
604 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
605
606 //
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);
614 fRx=fParam->GetPadRowRadii(fSector,row);
615
616
617 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
618 fZWidth = fParam->GetZWidth();
619 if (fSector < kNIS) {
620 fMaxPad = fParam->GetNPadsLow(row);
621 fSign = (fSector < kNIS/2) ? 1 : -1;
622 fPadLength = fParam->GetPadPitchLength(fSector,row);
623 fPadWidth = fParam->GetPadPitchWidth();
624 } else {
625 fMaxPad = fParam->GetNPadsUp(row);
626 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
627 fPadLength = fParam->GetPadPitchLength(fSector,row);
628 fPadWidth = fParam->GetPadPitchWidth();
629 }
630
631
632 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
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);
637
638 if (digarr.First()) //MI change
639 do {
640 Float_t dig=digarr.CurrentDigit();
641 if (dig<=fParam->GetZeroSup()) continue;
642 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
643 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
644 fBins[i*fMaxTime+j]=dig/gain;
645 } while (digarr.Next());
646 digarr.ExpandTrackBuffer();
647
648 FindClusters();
649
650 fOutput->Fill();
651 delete clrow;
652 nclusters+=fNcluster;
653 delete[] fBins;
654 delete[] fResBins;
655 }
656
657 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
658}
659
660void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
661{
662//-----------------------------------------------------------------
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
667//-----------------------------------------------------------------
668
669 if (!fOutput) {
670 Error("Digits2Clusters", "output tree not initialised");
671 return;
672 }
673
674 fRowDig = NULL;
675 AliTPCROC * roc = AliTPCROC::Instance();
676 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
677 AliTPCRawStream input(rawReader);
678 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
679 if (fEventHeader){
680 fTimeStamp = fEventHeader->Get("Timestamp");
681 fEventType = fEventHeader->Get("Type");
682 }
683
684
685 Int_t nclusters = 0;
686
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();
693 //
694 //alocate memory for sector - maximal case
695 //
696 Float_t** allBins = NULL;
697 Float_t** allBinsRes = NULL;
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 //
711 // Loop over sectors
712 //
713 for(fSector = 0; fSector < kNS; fSector++) {
714
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
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
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();
746 input.SetOldRCUFormat(fIsOldRCUFormat);
747 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
748 Int_t digCounter=0;
749 // Begin loop over altro data
750 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
751 Float_t gain =1;
752 Int_t lastPad=-1;
753 while (input.Next()) {
754 digCounter++;
755 if (input.GetSector() != fSector)
756 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
757
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));
763 //pad
764 Int_t iPad = input.GetPad();
765 if (iPad < 0 || iPad >= nPadsMax)
766 AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
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())
776 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
777 iTimeBin, 0, iTimeBin -1));
778 iTimeBin+=3;
779 //signal
780 Float_t signal = input.GetSignal();
781 if (!calcPedestal && signal <= zeroSup) continue;
782 allBins[iRow][iPad*fMaxTime+iTimeBin] = signal/gain;
783 allBins[iRow][iPad*fMaxTime+0] = 1; // pad with signal
784 } // End of the loop over altro data
785 //
786 //
787 // Now loop over rows and perform pedestal subtraction
788 if (digCounter==0) continue;
789 // if (fPedSubtraction) {
790 if (calcPedestal) {
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++) {
799 if (allBins[iRow][iPad*fMaxTime+0] !=1) continue; // no data
800 Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
801 //Float_t pedestal = TMath::Median(fMaxTime, p);
802 Int_t id[3] = {fSector, iRow, iPad-3};
803 Double_t rms=0;
804 Float_t pedestal = ProcesSignal(p, fMaxTime, id, rms);
805 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
806 allBins[iRow][iPad*fMaxTime+iTimeBin] -= pedestal;
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;
811 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
812 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
813 if (allBins[iRow][iPad*fMaxTime+iTimeBin] < 3.0*rms) // 3 sigma cut on RMS
814 allBins[iRow][iPad*fMaxTime+iTimeBin] = 0;
815 }
816 }
817 }
818 }
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);
829 fPadWidth = fParam->GetPadPitchWidth();
830 if (fSector < kNIS)
831 fMaxPad = fParam->GetNPadsLow(fRow);
832 else
833 fMaxPad = fParam->GetNPadsUp(fRow);
834 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
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
845 } // End of loop over sectors
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);
855
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 }
871 fBins[i+2*fMaxTime] = amp0;
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 }
881 fBins[(fMaxPad+3)*fMaxTime+i] = amp0;
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;
890 Float_t *b=&fBins[-1]+2*fMaxTime;
891 Int_t crtime = Int_t((fParam->GetZLength()-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
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 */
924}
925
926
927Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsOut){
928 //
929 // process signal on given pad - + streaming of additional information in special mode
930 //
931 // id[0] - sector
932 // id[1] - row
933 // id[2] - pad
934
935 //
936 // ESTIMATE pedestal and the noise
937 //
938 Int_t offset =100;
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;
945 //
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) {
951 max = signal[i];
952 maxPos = i;
953 }
954 if (signal[i]>kPedMax-1) continue;
955 histo[int(signal[i]+0.5)]++;
956 count0++;
957 }
958 //
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
964 //
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);
987 }
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));
1003 rmsOut = rms09;
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"<<
1012 "TimeStamp="<<fTimeStamp<<
1013 "EventType="<<fEventType<<
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])){
1039 TObjArray * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
1040 if (!sectorArray){
1041 Int_t npads =roc->GetNChannels(uid[0]);
1042 sectorArray = new TObjArray(npads);
1043 fAmplitudeHisto->AddAt(sectorArray, uid[0]);
1044 }
1045 Int_t position = uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
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();
1052 histo = new TH1F(hname, hname, 100, 5,100);
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++){
1058 histo->Fill(signal[i]);
1059 }
1060 }
1061 //
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 //
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
1079 "TimeStamp="<<fTimeStamp<<
1080 "EventType="<<fEventType<<
1081 "Sector="<<uid[0]<<
1082 "Row="<<uid[1]<<
1083 "Pad="<<uid[2]<<
1084 "Graph="<<graph<<
1085 "Max="<<max<<
1086 "MaxPos="<<maxPos<<
1087 //
1088 "Median="<<median<<
1089 "Mean="<<mean<<
1090 "RMS="<<rms<<
1091 "Mean06="<<mean06<<
1092 "RMS06="<<rms06<<
1093 "Mean09="<<mean09<<
1094 "RMS09="<<rms09<<
1095 "\n";
1096 if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1097 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1098 "TimeStamp="<<fTimeStamp<<
1099 "EventType="<<fEventType<<
1100 "Sector="<<uid[0]<<
1101 "Row="<<uid[1]<<
1102 "Pad="<<uid[2]<<
1103 "Graph="<<graph<<
1104 "Max="<<max<<
1105 "MaxPos="<<maxPos<<
1106 //
1107 "Median="<<median<<
1108 "Mean="<<mean<<
1109 "RMS="<<rms<<
1110 "Mean06="<<mean06<<
1111 "RMS06="<<rms06<<
1112 "Mean09="<<mean09<<
1113 "RMS09="<<rms09<<
1114 "\n";
1115 delete graph;
1116 }
1117
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;
1127 const Int_t kCemin=5; // range for the analysis of the ce signal +- channels from the peak
1128 const Int_t kCemax=5;
1129 for (Int_t i=nchannels-2; i>nchannels/2; i--){
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){
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;
1144 }
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;
1154 }
1155 }
1156 if (ceQmax&&ceQsum>ceSumThreshold) {
1157 ceTime/=ceQsum;
1158 (*fDebugStreamer)<<"Signalce"<<
1159 "TimeStamp="<<fTimeStamp<<
1160 "EventType="<<fEventType<<
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 }
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;
1183
1184 for (Int_t i=1; i<nchannels/4; i++){
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){
1187 ggmaxpos=i;
1188 if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
1189 break;
1190 }
1191 }
1192 if (ggmaxpos!=0){
1193 for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){
1194 if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
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"<<
1204 "TimeStamp="<<fTimeStamp<<
1205 "EventType="<<fEventType<<
1206 "Sector="<<uid[0]<<
1207 "Row="<<uid[1]<<
1208 "Pad="<<uid[2]<<
1209 "Max="<<ggQmax<<
1210 "Qsum="<<ggQsum<<
1211 "Time="<<ggTime<<
1212 "RMS06="<<rms06<<
1213 //
1214 "\n";
1215 }
1216 }
1217 // end of gg signal analysis
1218
1219
1220 delete [] dsignal;
1221 delete [] dtime;
1222 if (rms06>fRecoParam->GetMaxNoise()) return 1024+median; // sign noisy channel in debug mode
1223 return median;
1224}
1225
1226
1227
1228void AliTPCclustererMI::DumpHistos(){
1229 //
1230 // Dump histogram information
1231 //
1232 if (!fAmplitudeHisto) return;
1233 AliTPCROC * roc = AliTPCROC::Instance();
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;
1250 const UInt_t *indexes =roc->GetRowIndexes(isector);
1251 for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
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"<<
1260 "TimeStamp="<<fTimeStamp<<
1261 "EventType="<<fEventType<<
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}