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