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