]> 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
New trigger processing (towards an automatic procedure) and macro optimization (Cynthia)
[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// 1. The Input data for reconstruction - Options
22// 1.a Simulated data - TTree - invoked Digits2Clusters()
23// 1.b Raw data - Digits2Clusters(AliRawReader* rawReader);
24// 1.c HLT clusters - Digits2Clusters and Digits2Clusters(AliRawReader* rawReader)
25// invoke ReadHLTClusters()
26//
27// fUseHLTClusters - switches between different inputs
28// 1 -> only TPC raw/sim data
29// 2 -> if present TPC raw/sim data, otherwise HLT clusters
30// 3 -> only HLT clusters
31// 4 -> if present HLT clusters, otherwise TPC raw/sim data
32//
33// 2. The Output data
34// 2.a TTree with clusters - if SetOutput(TTree * tree) invoked
35// 2.b TObjArray - Faster option for HLT
36// 2.c TClonesArray - Faster option for HLT (smaller memory consumption), activate with fBClonesArray flag
37//
38// 3. Reconstruction setup
39// see AliTPCRecoParam for list of parameters
40// The reconstruction parameterization taken from the
41// AliTPCReconstructor::GetRecoParam()
42// Possible to setup it in reconstruction macro AliTPCReconstructor::SetRecoParam(...)
43//
44//
45//
46// Origin: Marian Ivanov
47//-------------------------------------------------------
48
49#include "Riostream.h"
50#include <TF1.h>
51#include <TFile.h>
52#include <TGraph.h>
53#include <TH1F.h>
54#include <TObjArray.h>
55#include <TClonesArray.h>
56#include <TRandom.h>
57#include <TTree.h>
58#include <TTreeStream.h>
59#include "TSystem.h"
60#include "TClass.h"
61
62#include "AliDigits.h"
63#include "AliLoader.h"
64#include "AliLog.h"
65#include "AliMathBase.h"
66#include "AliRawEventHeaderBase.h"
67#include "AliRawReader.h"
68#include "AliRunLoader.h"
69#include "AliSimDigits.h"
70#include "AliTPCCalPad.h"
71#include "AliTPCCalROC.h"
72#include "AliTPCClustersArray.h"
73#include "AliTPCClustersRow.h"
74#include "AliTPCParam.h"
75#include "AliTPCRawStream.h"
76#include "AliTPCRawStreamV3.h"
77#include "AliTPCRecoParam.h"
78#include "AliTPCReconstructor.h"
79#include "AliTPCcalibDB.h"
80#include "AliTPCclusterInfo.h"
81#include "AliTPCclusterMI.h"
82#include "AliTPCTransform.h"
83#include "AliTPCclustererMI.h"
84
85ClassImp(AliTPCclustererMI)
86
87
88
89AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
90 fBins(0),
91 fSigBins(0),
92 fNSigBins(0),
93 fLoop(0),
94 fMaxBin(0),
95 fMaxTime(1006), // 1000>940 so use 1000, add 3 virtual time bins before and 3 after
96 fMaxPad(0),
97 fSector(-1),
98 fRow(-1),
99 fSign(0),
100 fRx(0),
101 fPadWidth(0),
102 fPadLength(0),
103 fZWidth(0),
104 fPedSubtraction(kFALSE),
105 fEventHeader(0),
106 fTimeStamp(0),
107 fEventType(0),
108 fInput(0),
109 fOutput(0),
110 fOutputArray(0),
111 fOutputClonesArray(0),
112 fRowCl(0),
113 fRowDig(0),
114 fParam(0),
115 fNcluster(0),
116 fNclusters(0),
117 fDebugStreamer(0),
118 fRecoParam(0),
119 fBDumpSignal(kFALSE),
120 fBClonesArray(kFALSE),
121 fUseHLTClusters(1),
122 fAllBins(NULL),
123 fAllSigBins(NULL),
124 fAllNSigBins(NULL)
125{
126 //
127 // COSNTRUCTOR
128 // param - tpc parameters for given file
129 // recoparam - reconstruction parameters
130 //
131 fInput =0;
132 fParam = par;
133 if (recoParam) {
134 fRecoParam = recoParam;
135 }else{
136 //set default parameters if not specified
137 fRecoParam = AliTPCReconstructor::GetRecoParam();
138 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
139 }
140
141 if(AliTPCReconstructor::StreamLevel()>0) {
142 fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
143 }
144
145 // Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
146 fRowCl= new AliTPCClustersRow("AliTPCclusterMI");
147
148 // Non-persistent arrays
149 //
150 //alocate memory for sector - maximal case
151 //
152 AliTPCROC * roc = AliTPCROC::Instance();
153 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
154 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
155
156 fAllBins = new Float_t*[nRowsMax];
157 fAllSigBins = new Int_t*[nRowsMax];
158 fAllNSigBins = new Int_t[nRowsMax];
159 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
160 //
161 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
162 fAllBins[iRow] = new Float_t[maxBin];
163 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
164 fAllSigBins[iRow] = new Int_t[maxBin];
165 fAllNSigBins[iRow]=0;
166 }
167}
168//______________________________________________________________
169AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
170 :TObject(param),
171 fBins(0),
172 fSigBins(0),
173 fNSigBins(0),
174 fLoop(0),
175 fMaxBin(0),
176 fMaxTime(0),
177 fMaxPad(0),
178 fSector(-1),
179 fRow(-1),
180 fSign(0),
181 fRx(0),
182 fPadWidth(0),
183 fPadLength(0),
184 fZWidth(0),
185 fPedSubtraction(kFALSE),
186 fEventHeader(0),
187 fTimeStamp(0),
188 fEventType(0),
189 fInput(0),
190 fOutput(0),
191 fOutputArray(0),
192 fOutputClonesArray(0),
193 fRowCl(0),
194 fRowDig(0),
195 fParam(0),
196 fNcluster(0),
197 fNclusters(0),
198 fDebugStreamer(0),
199 fRecoParam(0),
200 fBDumpSignal(kFALSE),
201 fBClonesArray(kFALSE),
202 fUseHLTClusters(1),
203 fAllBins(NULL),
204 fAllSigBins(NULL),
205 fAllNSigBins(NULL)
206{
207 //
208 // dummy
209 //
210 fMaxBin = param.fMaxBin;
211}
212//______________________________________________________________
213AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
214{
215 //
216 // assignment operator - dummy
217 //
218 fMaxBin=param.fMaxBin;
219 return (*this);
220}
221//______________________________________________________________
222AliTPCclustererMI::~AliTPCclustererMI(){
223 //
224 //
225 //
226 if (fDebugStreamer) delete fDebugStreamer;
227 if (fOutputArray){
228 //fOutputArray->Delete();
229 delete fOutputArray;
230 }
231 if (fOutputClonesArray){
232 fOutputClonesArray->Delete();
233 delete fOutputClonesArray;
234 }
235
236 if (fRowCl) {
237 fRowCl->GetArray()->Delete();
238 delete fRowCl;
239 }
240
241 AliTPCROC * roc = AliTPCROC::Instance();
242 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
243 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
244 delete [] fAllBins[iRow];
245 delete [] fAllSigBins[iRow];
246 }
247 delete [] fAllBins;
248 delete [] fAllSigBins;
249 delete [] fAllNSigBins;
250}
251
252void AliTPCclustererMI::SetInput(TTree * tree)
253{
254 //
255 // set input tree with digits
256 //
257 fInput = tree;
258 if (!fInput->GetBranch("Segment")){
259 cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
260 fInput=0;
261 return;
262 }
263}
264
265void AliTPCclustererMI::SetOutput(TTree * tree)
266{
267 //
268 // Set the output tree
269 // If not set the ObjArray used - Option for HLT
270 //
271 if (!tree) return;
272 fOutput= tree;
273 AliTPCClustersRow clrow("AliTPCclusterMI");
274 AliTPCClustersRow *pclrow=&clrow;
275 fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
276}
277
278
279void AliTPCclustererMI::FillRow(){
280 //
281 // fill the output container -
282 // 2 Options possible
283 // Tree
284 // TObjArray
285 //
286 if (fOutput) fOutput->Fill();
287 if (!fOutput && !fBClonesArray){
288 //
289 if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
290 if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
291 }
292}
293
294Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
295 // sigma y2 = in digits - we don't know the angle
296 Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
297 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
298 (fPadWidth*fPadWidth);
299 Float_t sres = 0.25;
300 Float_t res = sd2+sres;
301 return res;
302}
303
304
305Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){
306 //sigma z2 = in digits - angle estimated supposing vertex constraint
307 Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
308 Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
309 Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
310 angular*=angular;
311 angular/=12.;
312 Float_t sres = fParam->GetZSigma()/fZWidth;
313 sres *=sres;
314 Float_t res = angular +sd2+sres;
315 return res;
316}
317
318void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
319AliTPCclusterMI &c)
320{
321 //
322 // k - Make cluster at position k
323 // bins - 2 D array of signals mapped to 1 dimensional array -
324 // max - the number of time bins er one dimension
325 // c - refernce to cluster to be filled
326 //
327 Int_t i0=k/max; //central pad
328 Int_t j0=k%max; //central time bin
329
330 // set pointers to data
331 //Int_t dummy[5] ={0,0,0,0,0};
332 Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
333 for (Int_t di=-2;di<=2;di++){
334 matrix[di+2] = &bins[k+di*max];
335 }
336 //build matrix with virtual charge
337 Float_t sigmay2= GetSigmaY2(j0);
338 Float_t sigmaz2= GetSigmaZ2(j0);
339
340 Float_t vmatrix[5][5];
341 vmatrix[2][2] = matrix[2][0];
342 c.SetType(0);
343 c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
344 for (Int_t di =-1;di <=1;di++)
345 for (Int_t dj =-1;dj <=1;dj++){
346 Float_t amp = matrix[di+2][dj];
347 if ( (amp<2) && (fLoop<2)){
348 // if under threshold - calculate virtual charge
349 Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
350 amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
351 if (amp>2) amp = 2;
352 vmatrix[2+di][2+dj]=amp;
353 vmatrix[2+2*di][2+2*dj]=0;
354 if ( (di*dj)!=0){
355 //DIAGONAL ELEMENTS
356 vmatrix[2+2*di][2+dj] =0;
357 vmatrix[2+di][2+2*dj] =0;
358 }
359 continue;
360 }
361 if (amp<4){
362 //if small amplitude - below 2 x threshold - don't consider other one
363 vmatrix[2+di][2+dj]=amp;
364 vmatrix[2+2*di][2+2*dj]=0; // don't take to the account next bin
365 if ( (di*dj)!=0){
366 //DIAGONAL ELEMENTS
367 vmatrix[2+2*di][2+dj] =0;
368 vmatrix[2+di][2+2*dj] =0;
369 }
370 continue;
371 }
372 //if bigger then take everything
373 vmatrix[2+di][2+dj]=amp;
374 vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
375 if ( (di*dj)!=0){
376 //DIAGONAL ELEMENTS
377 vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
378 vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
379 }
380 }
381
382
383
384 Float_t sumw=0;
385 Float_t sumiw=0;
386 Float_t sumi2w=0;
387 Float_t sumjw=0;
388 Float_t sumj2w=0;
389 //
390 for (Int_t i=-2;i<=2;i++)
391 for (Int_t j=-2;j<=2;j++){
392 Float_t amp = vmatrix[i+2][j+2];
393
394 sumw += amp;
395 sumiw += i*amp;
396 sumi2w += i*i*amp;
397 sumjw += j*amp;
398 sumj2w += j*j*amp;
399 }
400 //
401 Float_t meani = sumiw/sumw;
402 Float_t mi2 = sumi2w/sumw-meani*meani;
403 Float_t meanj = sumjw/sumw;
404 Float_t mj2 = sumj2w/sumw-meanj*meanj;
405 //
406 Float_t ry = mi2/sigmay2;
407 Float_t rz = mj2/sigmaz2;
408
409 //
410 if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
411 if ( ((ry <1.2) && (rz<1.2)) || (!fRecoParam->GetDoUnfold())) {
412 //
413 //if cluster looks like expected or Unfolding not switched on
414 //standard COG is used
415 //+1.2 deviation from expected sigma accepted
416 // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
417
418 meani +=i0;
419 meanj +=j0;
420 //set cluster parameters
421 c.SetQ(sumw);
422 c.SetPad(meani-2.5);
423 c.SetTimeBin(meanj-3);
424 c.SetSigmaY2(mi2);
425 c.SetSigmaZ2(mj2);
426 c.SetType(0);
427 AddCluster(c,(Float_t*)vmatrix,k);
428 return;
429 }
430 //
431 //unfolding when neccessary
432 //
433
434 Float_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3
435 Float_t dummy[7]={0,0,0,0,0,0};
436 for (Int_t di=-3;di<=3;di++){
437 matrix2[di+3] = &bins[k+di*max];
438 if ((k+di*max)<3) matrix2[di+3] = &dummy[3];
439 if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3];
440 }
441 Float_t vmatrix2[5][5];
442 Float_t sumu;
443 Float_t overlap;
444 UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
445 //
446 // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
447 meani +=i0;
448 meanj +=j0;
449 //set cluster parameters
450 c.SetQ(sumu);
451 c.SetPad(meani-2.5);
452 c.SetTimeBin(meanj-3);
453 c.SetSigmaY2(mi2);
454 c.SetSigmaZ2(mj2);
455 c.SetType(Char_t(overlap)+1);
456 AddCluster(c,(Float_t*)vmatrix,k);
457
458 //unfolding 2
459 meani-=i0;
460 meanj-=j0;
461}
462
463
464
465void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
466 Float_t & sumu, Float_t & overlap )
467{
468 //
469 //unfold cluster from input matrix
470 //data corresponding to cluster writen in recmatrix
471 //output meani and meanj
472
473 //take separatelly y and z
474
475 Float_t sum3i[7] = {0,0,0,0,0,0,0};
476 Float_t sum3j[7] = {0,0,0,0,0,0,0};
477
478 for (Int_t k =0;k<7;k++)
479 for (Int_t l = -1; l<=1;l++){
480 sum3i[k]+=matrix2[k][l];
481 sum3j[k]+=matrix2[l+3][k-3];
482 }
483 Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
484 //
485 //unfold y
486 Float_t sum3wi = 0; //charge minus overlap
487 Float_t sum3wio = 0; //full charge
488 Float_t sum3iw = 0; //sum for mean value
489 for (Int_t dk=-1;dk<=1;dk++){
490 sum3wio+=sum3i[dk+3];
491 if (dk==0){
492 sum3wi+=sum3i[dk+3];
493 }
494 else{
495 Float_t ratio =1;
496 if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
497 (sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 )){
498 Float_t xm2 = sum3i[-dk+3];
499 Float_t xm1 = sum3i[+3];
500 Float_t x1 = sum3i[2*dk+3];
501 Float_t x2 = sum3i[3*dk+3];
502 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
503 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
504 ratio = w11/(w11+w12);
505 for (Int_t dl=-1;dl<=1;dl++)
506 mratio[dk+1][dl+1] *= ratio;
507 }
508 Float_t amp = sum3i[dk+3]*ratio;
509 sum3wi+=amp;
510 sum3iw+= dk*amp;
511 }
512 }
513 meani = sum3iw/sum3wi;
514 Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
515
516
517
518 //unfold z
519 Float_t sum3wj = 0; //charge minus overlap
520 Float_t sum3wjo = 0; //full charge
521 Float_t sum3jw = 0; //sum for mean value
522 for (Int_t dk=-1;dk<=1;dk++){
523 sum3wjo+=sum3j[dk+3];
524 if (dk==0){
525 sum3wj+=sum3j[dk+3];
526 }
527 else{
528 Float_t ratio =1;
529 if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
530 (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
531 Float_t xm2 = sum3j[-dk+3];
532 Float_t xm1 = sum3j[+3];
533 Float_t x1 = sum3j[2*dk+3];
534 Float_t x2 = sum3j[3*dk+3];
535 Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
536 Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
537 ratio = w11/(w11+w12);
538 for (Int_t dl=-1;dl<=1;dl++)
539 mratio[dl+1][dk+1] *= ratio;
540 }
541 Float_t amp = sum3j[dk+3]*ratio;
542 sum3wj+=amp;
543 sum3jw+= dk*amp;
544 }
545 }
546 meanj = sum3jw/sum3wj;
547 Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
548 overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
549 sumu = (sum3wj+sum3wi)/2.;
550
551 if (overlap ==3) {
552 //if not overlap detected remove everything
553 for (Int_t di =-2; di<=2;di++)
554 for (Int_t dj =-2; dj<=2;dj++){
555 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
556 }
557 }
558 else{
559 for (Int_t di =-1; di<=1;di++)
560 for (Int_t dj =-1; dj<=1;dj++){
561 Float_t ratio =1;
562 if (mratio[di+1][dj+1]==1){
563 recmatrix[di+2][dj+2] = matrix2[3+di][dj];
564 if (TMath::Abs(di)+TMath::Abs(dj)>1){
565 recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
566 recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
567 }
568 recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
569 }
570 else
571 {
572 //if we have overlap in direction
573 recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
574 if (TMath::Abs(di)+TMath::Abs(dj)>1){
575 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
576 recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
577 //
578 ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
579 recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
580 }
581 else{
582 ratio = recmatrix[di+2][dj+2]/matrix2[3][0];
583 recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
584 }
585 }
586 }
587 }
588
589}
590
591Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
592{
593 //
594 // estimate max
595 Float_t sumteor= 0;
596 Float_t sumamp = 0;
597
598 for (Int_t di = -1;di<=1;di++)
599 for (Int_t dj = -1;dj<=1;dj++){
600 if (vmatrix[2+di][2+dj]>2){
601 Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
602 sumteor += teor*vmatrix[2+di][2+dj];
603 sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
604 }
605 }
606 Float_t max = sumamp/sumteor;
607 return max;
608}
609
610void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * /*matrix*/, Int_t /*pos*/){
611 //
612 //
613 // Transform cluster to the rotated global coordinata
614 // Assign labels to the cluster
615 // add the cluster to the array
616 // for more details - See AliTPCTranform::Transform(x,i,0,1)
617 Float_t meani = c.GetPad();
618 Float_t meanj = c.GetTimeBin();
619
620 Int_t ki = TMath::Nint(meani);
621 if (ki<0) ki=0;
622 if (ki>=fMaxPad) ki = fMaxPad-1;
623 Int_t kj = TMath::Nint(meanj);
624 if (kj<0) kj=0;
625 if (kj>=fMaxTime-3) kj=fMaxTime-4;
626 // ki and kj shifted as integers coordinata
627 if (fRowDig) {
628 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
629 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
630 c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
631 }
632
633 c.SetRow(fRow);
634 c.SetDetector(fSector);
635 Float_t s2 = c.GetSigmaY2();
636 Float_t w=fParam->GetPadPitchWidth(fSector);
637 c.SetSigmaY2(s2*w*w);
638 s2 = c.GetSigmaZ2();
639 c.SetSigmaZ2(s2*fZWidth*fZWidth);
640 //
641 //
642 //
643
644 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
645 if (!transform) {
646 AliFatal("Tranformations not in calibDB");
647 return;
648 }
649 transform->SetCurrentRecoParam((AliTPCRecoParam*)fRecoParam);
650 Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
651 Int_t i[1]={fSector};
652 transform->Transform(x,i,0,1);
653 c.SetX(x[0]);
654 c.SetY(x[1]);
655 c.SetZ(x[2]);
656 //
657 //
658 if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
659 c.SetType(-(c.GetType()+3)); //edge clusters
660 }
661 if (fLoop==2) c.SetType(100);
662 if (!AcceptCluster(&c)) return;
663
664 // select output
665 TClonesArray * arr = 0;
666 AliTPCclusterMI * cl = 0;
667
668 if(fBClonesArray==kFALSE) {
669 arr = fRowCl->GetArray();
670 cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
671 } else {
672 cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c);
673 }
674
675 // if (fRecoParam->DumpSignal() &&matrix ) {
676// Int_t nbins=0;
677// Float_t *graph =0;
678// if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
679// nbins = fMaxTime;
680// graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
681// }
682// AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
683// cl->SetInfo(info);
684// }
685 if (!fRecoParam->DumpSignal()) {
686 cl->SetInfo(0);
687 }
688
689 if (AliTPCReconstructor::StreamLevel()>1) {
690 Float_t xyz[3];
691 cl->GetGlobalXYZ(xyz);
692 (*fDebugStreamer)<<"Clusters"<<
693 "Cl.="<<cl<<
694 "gx="<<xyz[0]<<
695 "gy="<<xyz[1]<<
696 "gz="<<xyz[2]<<
697 "\n";
698 }
699
700 fNcluster++;
701}
702
703
704//_____________________________________________________________________________
705void AliTPCclustererMI::Digits2Clusters()
706{
707 //-----------------------------------------------------------------
708 // This is a simple cluster finder.
709 //-----------------------------------------------------------------
710
711 if (!fInput) {
712 Error("Digits2Clusters", "input tree not initialised");
713 return;
714 }
715 fRecoParam = AliTPCReconstructor::GetRecoParam();
716 if (!fRecoParam){
717 AliFatal("Can not get the reconstruction parameters");
718 }
719 if(AliTPCReconstructor::StreamLevel()>5) {
720 AliInfo("Parameter Dumps");
721 fParam->Dump();
722 fRecoParam->Dump();
723 }
724
725 //-----------------------------------------------------------------
726 // Use HLT clusters
727 //-----------------------------------------------------------------
728 fUseHLTClusters = fRecoParam->GetUseHLTClusters();
729
730 AliInfo(Form("Usage of HLT clusters in TPC reconstruction : %d",fUseHLTClusters));
731
732 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
733 AliInfo("Using HLT clusters for TPC off-line reconstruction");
734 fZWidth = fParam->GetZWidth();
735 Int_t iResult = ReadHLTClusters();
736
737 // HLT clusters present
738 if (iResult >= 0 && fNclusters > 0)
739 return;
740
741 // HLT clusters not present
742 if (iResult < 0 || fNclusters == 0) {
743 if (fUseHLTClusters == 3) {
744 AliError("No HLT clusters present, but requested.");
745 return;
746 }
747 else {
748 AliInfo("Now trying to read TPC RAW");
749 }
750 }
751 // Some other problem during cluster reading
752 else {
753 AliWarning("Some problem while unpacking of HLT clusters.");
754 return;
755 }
756 } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
757
758 //-----------------------------------------------------------------
759 // Run TPC off-line clusterer
760 //-----------------------------------------------------------------
761 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
762 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
763 AliSimDigits digarr, *dummy=&digarr;
764 fRowDig = dummy;
765 fInput->GetBranch("Segment")->SetAddress(&dummy);
766 Stat_t nentries = fInput->GetEntries();
767
768 fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after
769
770 Int_t nclusters = 0;
771
772 for (Int_t n=0; n<nentries; n++) {
773 fInput->GetEvent(n);
774 if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
775 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
776 continue;
777 }
778 Int_t row = fRow;
779 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
780 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
781 //
782
783 fRowCl->SetID(digarr.GetID());
784 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
785 fRx=fParam->GetPadRowRadii(fSector,row);
786
787
788 const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
789 fZWidth = fParam->GetZWidth();
790 if (fSector < kNIS) {
791 fMaxPad = fParam->GetNPadsLow(row);
792 fSign = (fSector < kNIS/2) ? 1 : -1;
793 fPadLength = fParam->GetPadPitchLength(fSector,row);
794 fPadWidth = fParam->GetPadPitchWidth();
795 } else {
796 fMaxPad = fParam->GetNPadsUp(row);
797 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
798 fPadLength = fParam->GetPadPitchLength(fSector,row);
799 fPadWidth = fParam->GetPadPitchWidth();
800 }
801
802
803 fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
804 fBins =new Float_t[fMaxBin];
805 fSigBins =new Int_t[fMaxBin];
806 fNSigBins = 0;
807 memset(fBins,0,sizeof(Float_t)*fMaxBin);
808
809 if (digarr.First()) //MI change
810 do {
811 Float_t dig=digarr.CurrentDigit();
812 if (dig<=fParam->GetZeroSup()) continue;
813 Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
814 Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
815 Int_t bin = i*fMaxTime+j;
816 if (gain>0){
817 fBins[bin]=dig/gain;
818 }else{
819 fBins[bin]=0;
820 }
821 fSigBins[fNSigBins++]=bin;
822 } while (digarr.Next());
823 digarr.ExpandTrackBuffer();
824
825 FindClusters(noiseROC);
826 FillRow();
827 fRowCl->GetArray()->Clear("C");
828 nclusters+=fNcluster;
829
830 delete[] fBins;
831 delete[] fSigBins;
832 }
833
834 Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
835
836 if (fUseHLTClusters == 2 && nclusters == 0) {
837 AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
838
839 fZWidth = fParam->GetZWidth();
840 ReadHLTClusters();
841 }
842}
843
844void AliTPCclustererMI::ProcessSectorData(){
845 //
846 // Process the data for the current sector
847 //
848
849 AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
850 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
851 AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector
852 AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector
853 //check the presence of the calibration
854 if (!noiseROC ||!pedestalROC ) {
855 AliError(Form("Missing calibration per sector\t%d\n",fSector));
856 return;
857 }
858 Int_t nRows=fParam->GetNRow(fSector);
859 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
860 Int_t zeroSup = fParam->GetZeroSup();
861 // if (calcPedestal) {
862 if (kFALSE ) {
863 for (Int_t iRow = 0; iRow < nRows; iRow++) {
864 Int_t maxPad = fParam->GetNPads(fSector, iRow);
865
866 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
867 //
868 // Temporary fix for data production - !!!! MARIAN
869 // The noise calibration should take mean and RMS - currently the Gaussian fit used
870 // In case of double peak - the pad should be rejected
871 //
872 // Line mean - if more than given digits over threshold - make a noise calculation
873 // and pedestal substration
874 if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue;
875 //
876 if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data
877 Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3];
878 //Float_t pedestal = TMath::Median(fMaxTime, p);
879 Int_t id[3] = {fSector, iRow, iPad-3};
880 // calib values
881 Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3);
882 Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
883 Double_t rmsEvent = rmsCalib;
884 Double_t pedestalEvent = pedestalCalib;
885 ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
886 if (rmsEvent<rmsCalib) rmsEvent = rmsCalib; // take worst scenario
887 if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
888
889 //
890 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
891 Int_t bin = iPad*fMaxTime+iTimeBin;
892 fAllBins[iRow][bin] -= pedestalEvent;
893 if (iTimeBin < fRecoParam->GetFirstBin())
894 fAllBins[iRow][bin] = 0;
895 if (iTimeBin > fRecoParam->GetLastBin())
896 fAllBins[iRow][bin] = 0;
897 if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
898 fAllBins[iRow][bin] = 0;
899 if (fAllBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS
900 fAllBins[iRow][bin] = 0;
901 if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
902 }
903 }
904 }
905 }
906
907 if (AliTPCReconstructor::StreamLevel()>5) {
908 for (Int_t iRow = 0; iRow < nRows; iRow++) {
909 Int_t maxPad = fParam->GetNPads(fSector,iRow);
910
911 for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
912 for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
913 Int_t bin = iPad*fMaxTime+iTimeBin;
914 Float_t signal = fAllBins[iRow][bin];
915 if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
916 Double_t x[]={iRow,iPad-3,iTimeBin-3};
917 Int_t i[]={fSector};
918 AliTPCTransform trafo;
919 trafo.Transform(x,i,0,1);
920 Double_t gx[3]={x[0],x[1],x[2]};
921 trafo.RotatedGlobal2Global(fSector,gx);
922 // fAllSigBins[iRow][fAllNSigBins[iRow]++]
923 Int_t rowsigBins = fAllNSigBins[iRow];
924 Int_t first=fAllSigBins[iRow][0];
925 Int_t last= 0;
926 // if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1];
927
928 if (AliTPCReconstructor::StreamLevel()>5) {
929 (*fDebugStreamer)<<"Digits"<<
930 "sec="<<fSector<<
931 "row="<<iRow<<
932 "pad="<<iPad<<
933 "time="<<iTimeBin<<
934 "sig="<<signal<<
935 "x="<<x[0]<<
936 "y="<<x[1]<<
937 "z="<<x[2]<<
938 "gx="<<gx[0]<<
939 "gy="<<gx[1]<<
940 "gz="<<gx[2]<<
941 //
942 "rowsigBins="<<rowsigBins<<
943 "first="<<first<<
944 "last="<<last<<
945 "\n";
946 }
947 }
948 }
949 }
950 }
951 }
952
953 // Now loop over rows and find clusters
954 for (fRow = 0; fRow < nRows; fRow++) {
955 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
956 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
957
958 fRx = fParam->GetPadRowRadii(fSector, fRow);
959 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
960 fPadWidth = fParam->GetPadPitchWidth();
961 fMaxPad = fParam->GetNPads(fSector,fRow);
962 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
963
964 fBins = fAllBins[fRow];
965 fSigBins = fAllSigBins[fRow];
966 fNSigBins = fAllNSigBins[fRow];
967
968 FindClusters(noiseROC);
969
970 FillRow();
971 if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C");
972 fNclusters += fNcluster;
973
974 } // End of loop to find clusters
975}
976
977
978void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
979{
980//-----------------------------------------------------------------
981// This is a cluster finder for the TPC raw data.
982// The method assumes NO ordering of the altro channels.
983// The pedestal subtraction can be switched on and off
984// using an option of the TPC reconstructor
985//-----------------------------------------------------------------
986 fRecoParam = AliTPCReconstructor::GetRecoParam();
987 if (!fRecoParam){
988 AliFatal("Can not get the reconstruction parameters");
989 }
990 if(AliTPCReconstructor::StreamLevel()>5) {
991 AliInfo("Parameter Dumps");
992 fParam->Dump();
993 fRecoParam->Dump();
994 }
995 fRowDig = NULL;
996
997 //-----------------------------------------------------------------
998 // Use HLT clusters
999 //-----------------------------------------------------------------
1000 fUseHLTClusters = fRecoParam->GetUseHLTClusters();
1001
1002 AliInfo(Form("Usage of HLT clusters in TPC reconstruction : %d",fUseHLTClusters));
1003
1004 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
1005 AliInfo("Using HLT clusters for TPC off-line reconstruction");
1006 fZWidth = fParam->GetZWidth();
1007 Int_t iResult = ReadHLTClusters();
1008
1009 // HLT clusters present
1010 if (iResult >= 0 && fNclusters > 0)
1011 return;
1012
1013 // HLT clusters not present
1014 if (iResult < 0 || fNclusters == 0) {
1015 if (fUseHLTClusters == 3) {
1016 AliError("No HLT clusters present, but requested.");
1017 return;
1018 }
1019 else {
1020 AliInfo("Now trying to read TPC RAW");
1021 }
1022 }
1023 // Some other problem during cluster reading
1024 else {
1025 AliWarning("Some problem while unpacking of HLT clusters.");
1026 return;
1027 }
1028 } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
1029
1030 //-----------------------------------------------------------------
1031 // Run TPC off-line clusterer
1032 //-----------------------------------------------------------------
1033 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
1034 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1035 //
1036 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
1037 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1038 if (fEventHeader){
1039 fTimeStamp = fEventHeader->Get("Timestamp");
1040 fEventType = fEventHeader->Get("Type");
1041 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1042 transform->SetCurrentTimeStamp(fTimeStamp);
1043 transform->SetCurrentRun(rawReader->GetRunNumber());
1044 }
1045
1046 // creaate one TClonesArray for all clusters
1047 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
1048 // reset counter
1049 fNclusters = 0;
1050
1051 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
1052// const Int_t kNIS = fParam->GetNInnerSector();
1053// const Int_t kNOS = fParam->GetNOuterSector();
1054// const Int_t kNS = kNIS + kNOS;
1055 fZWidth = fParam->GetZWidth();
1056 Int_t zeroSup = fParam->GetZeroSup();
1057 //
1058 // Clean-up
1059 //
1060 AliTPCROC * roc = AliTPCROC::Instance();
1061 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1062 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1063 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1064 //
1065 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
1066 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1067 fAllNSigBins[iRow]=0;
1068 }
1069
1070 Int_t prevSector=-1;
1071 rawReader->Reset();
1072 Int_t digCounter=0;
1073 //
1074 // Loop over DDLs
1075 //
1076 const Int_t kNIS = fParam->GetNInnerSector();
1077 const Int_t kNOS = fParam->GetNOuterSector();
1078 const Int_t kNS = kNIS + kNOS;
1079
1080 for(fSector = 0; fSector < kNS; fSector++) {
1081
1082 Int_t nRows = 0;
1083 Int_t nDDLs = 0, indexDDL = 0;
1084 if (fSector < kNIS) {
1085 nRows = fParam->GetNRowLow();
1086 fSign = (fSector < kNIS/2) ? 1 : -1;
1087 nDDLs = 2;
1088 indexDDL = fSector * 2;
1089 }
1090 else {
1091 nRows = fParam->GetNRowUp();
1092 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1093 nDDLs = 4;
1094 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1095 }
1096
1097 // load the raw data for corresponding DDLs
1098 rawReader->Reset();
1099 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1100
1101 while (input.NextDDL()){
1102 if (input.GetSector() != fSector)
1103 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1104
1105 //Int_t nRows = fParam->GetNRow(fSector);
1106
1107 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
1108 // Begin loop over altro data
1109 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1110 Float_t gain =1;
1111
1112 //loop over pads
1113 while ( input.NextChannel() ) {
1114 Int_t iRow = input.GetRow();
1115 if (iRow < 0){
1116 continue;
1117 }
1118 if (iRow >= nRows){
1119 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1120 iRow, 0, nRows -1));
1121 continue;
1122 }
1123 //pad
1124 Int_t iPad = input.GetPad();
1125 if (iPad < 0 || iPad >= nPadsMax) {
1126 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1127 iPad, 0, nPadsMax-1));
1128 continue;
1129 }
1130 gain = gainROC->GetValue(iRow,iPad);
1131 iPad+=3;
1132
1133 //loop over bunches
1134 while ( input.NextBunch() ){
1135 Int_t startTbin = (Int_t)input.GetStartTimeBin();
1136 Int_t bunchlength = (Int_t)input.GetBunchLength();
1137 const UShort_t *sig = input.GetSignals();
1138 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1139 Int_t iTimeBin=startTbin-iTime;
1140 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1141 continue;
1142 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1143 iTimeBin, 0, iTimeBin -1));
1144 }
1145 iTimeBin+=3;
1146 //signal
1147 Float_t signal=(Float_t)sig[iTime];
1148 if (!calcPedestal && signal <= zeroSup) continue;
1149
1150 if (!calcPedestal) {
1151 Int_t bin = iPad*fMaxTime+iTimeBin;
1152 if (gain>0){
1153 fAllBins[iRow][bin] = signal/gain;
1154 }else{
1155 fAllBins[iRow][bin] =0;
1156 }
1157 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
1158 }else{
1159 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1160 }
1161 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
1162
1163 // Temporary
1164 digCounter++;
1165 }// end loop signals in bunch
1166 }// end loop bunches
1167 } // end loop pads
1168 //
1169 //
1170 //
1171 //
1172 // Now loop over rows and perform pedestal subtraction
1173 if (digCounter==0) continue;
1174 } // End of loop over sectors
1175 //process last sector
1176 if ( digCounter>0 ){
1177 ProcessSectorData();
1178 for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1179 Int_t maxPad = fParam->GetNPads(fSector,iRow);
1180 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
1181 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1182 fAllNSigBins[iRow] = 0;
1183 }
1184 prevSector=fSector;
1185 digCounter=0;
1186 }
1187 }
1188
1189 if (rawReader->GetEventId() && fOutput ){
1190 Info("Digits2Clusters", "File %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1191 }
1192
1193 if(rawReader->GetEventId()) {
1194 Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1195 }
1196
1197 if(fBClonesArray) {
1198 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1199 }
1200
1201 if (fUseHLTClusters == 2 && fNclusters == 0) {
1202 AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
1203
1204 fZWidth = fParam->GetZWidth();
1205 ReadHLTClusters();
1206 }
1207}
1208
1209
1210
1211
1212
1213void AliTPCclustererMI::Digits2ClustersOld
1214(AliRawReader* rawReader)
1215{
1216//-----------------------------------------------------------------
1217// This is a cluster finder for the TPC raw data.
1218// The method assumes NO ordering of the altro channels.
1219// The pedestal subtraction can be switched on and off
1220// using an option of the TPC reconstructor
1221//-----------------------------------------------------------------
1222 fRecoParam = AliTPCReconstructor::GetRecoParam();
1223 if (!fRecoParam){
1224 AliFatal("Can not get the reconstruction parameters");
1225 }
1226 if(AliTPCReconstructor::StreamLevel()>5) {
1227 AliInfo("Parameter Dumps");
1228 fParam->Dump();
1229 fRecoParam->Dump();
1230 }
1231 fRowDig = NULL;
1232
1233 AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
1234 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1235 //
1236 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
1237 fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1238 if (fEventHeader){
1239 fTimeStamp = fEventHeader->Get("Timestamp");
1240 fEventType = fEventHeader->Get("Type");
1241 }
1242
1243 // creaate one TClonesArray for all clusters
1244 if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
1245 // reset counter
1246 fNclusters = 0;
1247
1248 fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
1249 const Int_t kNIS = fParam->GetNInnerSector();
1250 const Int_t kNOS = fParam->GetNOuterSector();
1251 const Int_t kNS = kNIS + kNOS;
1252 fZWidth = fParam->GetZWidth();
1253 Int_t zeroSup = fParam->GetZeroSup();
1254 //
1255 // Clean-up
1256 //
1257
1258 AliTPCROC * roc = AliTPCROC::Instance();
1259 Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1260 Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1261 for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1262 //
1263 Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after
1264 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1265 fAllNSigBins[iRow]=0;
1266 }
1267 //
1268 // Loop over sectors
1269 //
1270 for(fSector = 0; fSector < kNS; fSector++) {
1271
1272 Int_t nRows = 0;
1273 Int_t nDDLs = 0, indexDDL = 0;
1274 if (fSector < kNIS) {
1275 nRows = fParam->GetNRowLow();
1276 fSign = (fSector < kNIS/2) ? 1 : -1;
1277 nDDLs = 2;
1278 indexDDL = fSector * 2;
1279 }
1280 else {
1281 nRows = fParam->GetNRowUp();
1282 fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1283 nDDLs = 4;
1284 indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1285 }
1286
1287 // load the raw data for corresponding DDLs
1288 rawReader->Reset();
1289 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1290
1291 // select only good sector
1292 if (!input.Next()) continue;
1293 if(input.GetSector() != fSector) continue;
1294
1295 AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector
1296
1297 for (Int_t iRow = 0; iRow < nRows; iRow++) {
1298 Int_t maxPad;
1299 if (fSector < kNIS)
1300 maxPad = fParam->GetNPadsLow(iRow);
1301 else
1302 maxPad = fParam->GetNPadsUp(iRow);
1303
1304 Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after
1305 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1306 fAllNSigBins[iRow] = 0;
1307 }
1308
1309 Int_t digCounter=0;
1310 // Begin loop over altro data
1311 Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1312 Float_t gain =1;
1313 Int_t lastPad=-1;
1314
1315 input.Reset();
1316 while (input.Next()) {
1317 if (input.GetSector() != fSector)
1318 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1319
1320
1321 Int_t iRow = input.GetRow();
1322 if (iRow < 0){
1323 continue;
1324 }
1325
1326 if (iRow < 0 || iRow >= nRows){
1327 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1328 iRow, 0, nRows -1));
1329 continue;
1330 }
1331 //pad
1332 Int_t iPad = input.GetPad();
1333 if (iPad < 0 || iPad >= nPadsMax) {
1334 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1335 iPad, 0, nPadsMax-1));
1336 continue;
1337 }
1338 if (iPad!=lastPad){
1339 gain = gainROC->GetValue(iRow,iPad);
1340 lastPad = iPad;
1341 }
1342 iPad+=3;
1343 //time
1344 Int_t iTimeBin = input.GetTime();
1345 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1346 continue;
1347 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1348 iTimeBin, 0, iTimeBin -1));
1349 }
1350 iTimeBin+=3;
1351
1352 //signal
1353 Float_t signal = input.GetSignal();
1354 if (!calcPedestal && signal <= zeroSup) continue;
1355
1356 if (!calcPedestal) {
1357 Int_t bin = iPad*fMaxTime+iTimeBin;
1358 if (gain>0){
1359 fAllBins[iRow][bin] = signal/gain;
1360 }else{
1361 fAllBins[iRow][bin] =0;
1362 }
1363 fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
1364 }else{
1365 fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1366 }
1367 fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal
1368
1369 // Temporary
1370 digCounter++;
1371 } // End of the loop over altro data
1372 //
1373 //
1374 //
1375 //
1376 // Now loop over rows and perform pedestal subtraction
1377 if (digCounter==0) continue;
1378 ProcessSectorData();
1379 } // End of loop over sectors
1380
1381 if (rawReader->GetEventId() && fOutput ){
1382 Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1383 }
1384
1385 if(rawReader->GetEventId()) {
1386 Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1387 }
1388
1389 if(fBClonesArray) {
1390 //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1391 }
1392}
1393
1394void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
1395{
1396
1397 //
1398 // add virtual charge at the edge
1399 //
1400 Double_t kMaxDumpSize = 500000;
1401 if (!fOutput) {
1402 fBDumpSignal =kFALSE;
1403 }else{
1404 if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag
1405 }
1406
1407 fNcluster=0;
1408 fLoop=1;
1409 Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
1410 Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs();
1411 Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1412 Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs();
1413 Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma();
1414 Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1415 Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma();
1416 Int_t useOnePadCluster = fRecoParam->GetUseOnePadCluster();
1417 for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1418 Int_t i = fSigBins[iSig];
1419 if (i%fMaxTime<=crtime) continue;
1420 Float_t *b = &fBins[i];
1421 //absolute custs
1422 if (b[0]<minMaxCutAbs) continue; //threshold for maxima
1423 //
1424 if (useOnePadCluster==0){
1425 if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1426 if (b[-1]+b[1]<=0) continue; // cut on isolated clusters
1427 if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1428 }
1429 //
1430 if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue; //threshold for up down (TRF)
1431 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue; //threshold for left right (PRF)
1432 if (!IsMaximum(*b,fMaxTime,b)) continue;
1433 //
1434 Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
1435 if (noise>fRecoParam->GetMaxNoise()) continue;
1436 // sigma cuts
1437 if (b[0]<minMaxCutSigma*noise) continue; //threshold form maxima
1438 if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue; //threshold for up town TRF
1439 if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue; //threshold for left right (PRF)
1440
1441 AliTPCclusterMI c; // default cosntruction without info
1442 Int_t dummy=0;
1443 MakeCluster(i, fMaxTime, fBins, dummy,c);
1444
1445 //}
1446 }
1447}
1448
1449Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1450 //
1451 // Currently hack to filter digital noise (15.06.2008)
1452 // To be parameterized in the AliTPCrecoParam
1453 // More inteligent way to be used in future
1454 // Acces to the proper pedestal file needed
1455 //
1456 if (cl->GetMax()<400) return kTRUE;
1457 Double_t ratio = cl->GetQ()/cl->GetMax();
1458 if (cl->GetMax()>700){
1459 if ((ratio - int(ratio)>0.8)) return kFALSE;
1460 }
1461 if ((ratio - int(ratio)<0.95)) return kTRUE;
1462 return kFALSE;
1463}
1464
1465
1466Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
1467 //
1468 // process signal on given pad - + streaming of additional information in special mode
1469 //
1470 // id[0] - sector
1471 // id[1] - row
1472 // id[2] - pad
1473
1474 //
1475 // ESTIMATE pedestal and the noise
1476 //
1477 const Int_t kPedMax = 100;
1478 Float_t max = 0;
1479 Float_t maxPos = 0;
1480 Int_t median = -1;
1481 Int_t count0 = 0;
1482 Int_t count1 = 0;
1483 Float_t rmsCalib = rmsEvent; // backup initial value ( from calib)
1484 Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1485 Int_t firstBin = fRecoParam->GetFirstBin();
1486 //
1487 UShort_t histo[kPedMax];
1488 //memset(histo,0,kPedMax*sizeof(UShort_t));
1489 for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
1490 for (Int_t i=0; i<fMaxTime; i++){
1491 if (signal[i]<=0) continue;
1492 if (signal[i]>max && i>firstBin) {
1493 max = signal[i];
1494 maxPos = i;
1495 }
1496 if (signal[i]>kPedMax-1) continue;
1497 histo[int(signal[i]+0.5)]++;
1498 count0++;
1499 }
1500 //
1501 for (Int_t i=1; i<kPedMax; i++){
1502 if (count1<count0*0.5) median=i;
1503 count1+=histo[i];
1504 }
1505 // truncated mean
1506 //
1507 Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
1508 Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ;
1509 Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ;
1510 //
1511 for (Int_t idelta=1; idelta<10; idelta++){
1512 if (median-idelta<=0) continue;
1513 if (median+idelta>kPedMax) continue;
1514 if (count06<0.6*count1){
1515 count06+=histo[median-idelta];
1516 mean06 +=histo[median-idelta]*(median-idelta);
1517 rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1518 count06+=histo[median+idelta];
1519 mean06 +=histo[median+idelta]*(median+idelta);
1520 rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1521 }
1522 if (count09<0.9*count1){
1523 count09+=histo[median-idelta];
1524 mean09 +=histo[median-idelta]*(median-idelta);
1525 rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta);
1526 count09+=histo[median+idelta];
1527 mean09 +=histo[median+idelta]*(median+idelta);
1528 rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta);
1529 }
1530 if (count10<0.95*count1){
1531 count10+=histo[median-idelta];
1532 mean +=histo[median-idelta]*(median-idelta);
1533 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
1534 count10+=histo[median+idelta];
1535 mean +=histo[median+idelta]*(median+idelta);
1536 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
1537 }
1538 }
1539 if (count10) {
1540 mean /=count10;
1541 rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1542 }
1543 if (count06) {
1544 mean06/=count06;
1545 rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1546 }
1547 if (count09) {
1548 mean09/=count09;
1549 rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1550 }
1551 rmsEvent = rms09;
1552 //
1553 pedestalEvent = median;
1554 if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1555 //
1556 UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1557 //
1558 // Dump mean signal info
1559 //
1560 if (AliTPCReconstructor::StreamLevel()>0) {
1561 (*fDebugStreamer)<<"Signal"<<
1562 "TimeStamp="<<fTimeStamp<<
1563 "EventType="<<fEventType<<
1564 "Sector="<<uid[0]<<
1565 "Row="<<uid[1]<<
1566 "Pad="<<uid[2]<<
1567 "Max="<<max<<
1568 "MaxPos="<<maxPos<<
1569 //
1570 "Median="<<median<<
1571 "Mean="<<mean<<
1572 "RMS="<<rms<<
1573 "Mean06="<<mean06<<
1574 "RMS06="<<rms06<<
1575 "Mean09="<<mean09<<
1576 "RMS09="<<rms09<<
1577 "RMSCalib="<<rmsCalib<<
1578 "PedCalib="<<pedestalCalib<<
1579 "\n";
1580 }
1581 //
1582 // fill pedestal histogram
1583 //
1584 //
1585 //
1586 //
1587 Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped
1588 Float_t *dsignal = new Float_t[nchannels];
1589 Float_t *dtime = new Float_t[nchannels];
1590 for (Int_t i=0; i<nchannels; i++){
1591 dtime[i] = i;
1592 dsignal[i] = signal[i];
1593 }
1594
1595 //
1596 // Big signals dumping
1597 //
1598 if (AliTPCReconstructor::StreamLevel()>0) {
1599 if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin())
1600 (*fDebugStreamer)<<"SignalB"<< // pads with signal
1601 "TimeStamp="<<fTimeStamp<<
1602 "EventType="<<fEventType<<
1603 "Sector="<<uid[0]<<
1604 "Row="<<uid[1]<<
1605 "Pad="<<uid[2]<<
1606 // "Graph="<<graph<<
1607 "Max="<<max<<
1608 "MaxPos="<<maxPos<<
1609 //
1610 "Median="<<median<<
1611 "Mean="<<mean<<
1612 "RMS="<<rms<<
1613 "Mean06="<<mean06<<
1614 "RMS06="<<rms06<<
1615 "Mean09="<<mean09<<
1616 "RMS09="<<rms09<<
1617 "\n";
1618 }
1619
1620 delete [] dsignal;
1621 delete [] dtime;
1622 if (rms06>fRecoParam->GetMaxNoise()) {
1623 pedestalEvent+=1024.;
1624 return 1024+median; // sign noisy channel in debug mode
1625 }
1626 return median;
1627}
1628
1629Int_t AliTPCclustererMI::ReadHLTClusters()
1630{
1631 //
1632 // read HLT clusters instead of off line custers,
1633 // used in Digits2Clusters
1634 //
1635
1636 TObject* pClusterAccess=NULL;
1637 TClass* pCl=NULL;
1638 ROOT::NewFunc_t pNewFunc=NULL;
1639 do {
1640 pCl=TClass::GetClass("AliHLTTPCClusterAccessHLTOUT");
1641 } while (!pCl && gSystem->Load("libAliHLTTPC.so")==0);
1642 if (!pCl || (pNewFunc=pCl->GetNew())==NULL) {
1643 AliError("can not load class description of AliHLTTPCClusterAccessHLTOUT, aborting ...");
1644 return -1;
1645 }
1646
1647 void* p=(*pNewFunc)(NULL);
1648 if (!p) {
1649 AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT");
1650 return -2;
1651 }
1652 pClusterAccess=reinterpret_cast<TObject*>(p);
1653 if (!pClusterAccess) {
1654 AliError("instance not of type TObject");
1655 return -3 ;
1656 }
1657
1658 const Int_t kNIS = fParam->GetNInnerSector();
1659 const Int_t kNOS = fParam->GetNOuterSector();
1660 const Int_t kNS = kNIS + kNOS;
1661 fNclusters = 0;
1662
1663 for(fSector = 0; fSector < kNS; fSector++) {
1664
1665 Int_t iResult = 1;
1666 TString param("sector="); param+=fSector;
1667 pClusterAccess->Clear();
1668 pClusterAccess->Execute("read", param, &iResult);
1669 if (iResult < 0) {
1670 return iResult;
1671 AliError("HLT Clusters can not be found");
1672 }
1673
1674 if (pClusterAccess->FindObject("clusterarray")==NULL) {
1675 AliError("HLT clusters requested, but not cluster array not present");
1676 return -4;
1677 }
1678
1679 TClonesArray* clusterArray=dynamic_cast<TClonesArray*>(pClusterAccess->FindObject("clusterarray"));
1680 if (!clusterArray) {
1681 AliError("HLT cluster array is not of class type TClonesArray");
1682 return -5;
1683 }
1684
1685 AliDebug(4,Form("Reading %d clusters from HLT for sector %d", clusterArray->GetEntriesFast(), fSector));
1686
1687 Int_t nClusterSector=0;
1688 Int_t nRows=fParam->GetNRow(fSector);
1689
1690 for (fRow = 0; fRow < nRows; fRow++) {
1691 fRowCl->SetID(fParam->GetIndex(fSector, fRow));
1692 if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
1693 fNcluster=0; // reset clusters per row
1694
1695 fRx = fParam->GetPadRowRadii(fSector, fRow);
1696 fPadLength = fParam->GetPadPitchLength(fSector, fRow);
1697 fPadWidth = fParam->GetPadPitchWidth();
1698 fMaxPad = fParam->GetNPads(fSector,fRow);
1699 fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
1700
1701 fBins = fAllBins[fRow];
1702 fSigBins = fAllSigBins[fRow];
1703 fNSigBins = fAllNSigBins[fRow];
1704
1705 for (Int_t i=0; i<clusterArray->GetEntriesFast(); i++) {
1706 if (!clusterArray->At(i))
1707 continue;
1708
1709 AliTPCclusterMI* cluster=dynamic_cast<AliTPCclusterMI*>(clusterArray->At(i));
1710 if (!cluster) continue;
1711 if (cluster->GetRow()!=fRow) continue;
1712 nClusterSector++;
1713 AddCluster(*cluster, NULL, 0);
1714 }
1715
1716 FillRow();
1717 fRowCl->GetArray()->Clear("c");
1718
1719 } // for (fRow = 0; fRow < nRows; fRow++) {
1720 if (nClusterSector!=clusterArray->GetEntriesFast()) {
1721 AliError(Form("Failed to read %d out of %d HLT clusters",
1722 clusterArray->GetEntriesFast()-nClusterSector,
1723 clusterArray->GetEntriesFast()));
1724 }
1725 fNclusters+=nClusterSector;
1726 } // for(fSector = 0; fSector < kNS; fSector++) {
1727
1728 delete pClusterAccess;
1729
1730 Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters);
1731
1732 return 0;
1733}