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