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