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