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