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