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