16145915a6ad71dd09292498f0835945c873fdcd
[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 //   Origin: Marian Ivanov 
22 //-------------------------------------------------------
23
24 #include "Riostream.h"
25 #include <TF1.h>
26 #include <TFile.h>
27 #include <TGraph.h>
28 #include <TH1F.h>
29 #include <TObjArray.h>
30 #include <TRandom.h>
31 #include <TTree.h>
32 #include <TTreeStream.h>
33
34 #include "AliDigits.h"
35 #include "AliLoader.h"
36 #include "AliLog.h"
37 #include "AliMathBase.h"
38 #include "AliRawEventHeaderBase.h"
39 #include "AliRawReader.h"
40 #include "AliRunLoader.h"
41 #include "AliSimDigits.h"
42 #include "AliTPCCalPad.h"
43 #include "AliTPCCalROC.h"
44 #include "AliTPCClustersArray.h"
45 #include "AliTPCClustersRow.h"
46 #include "AliTPCParam.h"
47 #include "AliTPCRawStream.h"
48 #include "AliTPCRecoParam.h"
49 #include "AliTPCReconstructor.h"
50 #include "AliTPCcalibDB.h"
51 #include "AliTPCclusterInfo.h"
52 #include "AliTPCclusterMI.h"
53 #include "AliTPCTransform.h"
54 #include "AliTPCclustererMI.h"
55
56 ClassImp(AliTPCclustererMI)
57
58
59
60 AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
61   fBins(0),
62   fSigBins(0),
63   fNSigBins(0),
64   fLoop(0),
65   fMaxBin(0),
66   fMaxTime(0),
67   fMaxPad(0),
68   fSector(-1),
69   fRow(-1),
70   fSign(0),
71   fRx(0),
72   fPadWidth(0),
73   fPadLength(0),
74   fZWidth(0),
75   fPedSubtraction(kFALSE),
76   fIsOldRCUFormat(kFALSE),
77   fEventHeader(0),
78   fTimeStamp(0),
79   fEventType(0),
80   fInput(0),
81   fOutput(0),
82   fRowCl(0),
83   fRowDig(0),
84   fParam(0),
85   fNcluster(0),
86   fDebugStreamer(0),
87   fRecoParam(0),
88   fBDumpSignal(kFALSE)
89 {
90   //
91   // COSNTRUCTOR
92   // param     - tpc parameters for given file
93   // recoparam - reconstruction parameters 
94   //
95   fIsOldRCUFormat = kFALSE;
96   fInput =0;
97   fOutput=0;
98   fParam = par;
99   if (recoParam) {
100     fRecoParam = recoParam;
101   }else{
102     //set default parameters if not specified
103     fRecoParam = AliTPCReconstructor::GetRecoParam();
104     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
105   }
106   fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
107   Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
108 }
109 //______________________________________________________________
110 AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
111               :TObject(param),
112   fBins(0),
113   fSigBins(0),
114   fNSigBins(0),
115   fLoop(0),
116   fMaxBin(0),
117   fMaxTime(0),
118   fMaxPad(0),
119   fSector(-1),
120   fRow(-1),
121   fSign(0),
122   fRx(0),
123   fPadWidth(0),
124   fPadLength(0),
125   fZWidth(0),
126   fPedSubtraction(kFALSE),
127   fIsOldRCUFormat(kFALSE),
128   fEventHeader(0),
129   fTimeStamp(0),
130   fEventType(0),
131   fInput(0),
132   fOutput(0),
133   fRowCl(0),
134   fRowDig(0),
135   fParam(0),
136   fNcluster(0),
137   fDebugStreamer(0),
138   fRecoParam(0),
139   fBDumpSignal(kFALSE)
140 {
141   //
142   // dummy
143   //
144   fMaxBin = param.fMaxBin;
145 }
146 //______________________________________________________________
147 AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
148 {
149   //
150   // assignment operator - dummy
151   //
152   fMaxBin=param.fMaxBin;
153   return (*this);
154 }
155 //______________________________________________________________
156 AliTPCclustererMI::~AliTPCclustererMI(){
157   if (fDebugStreamer) delete fDebugStreamer;
158 }
159
160 void AliTPCclustererMI::SetInput(TTree * tree)
161 {
162   //
163   // set input tree with digits
164   //
165   fInput = tree;  
166   if  (!fInput->GetBranch("Segment")){
167     cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
168     fInput=0;
169     return;
170   }
171 }
172
173 void AliTPCclustererMI::SetOutput(TTree * tree) 
174 {
175   //
176   //
177   fOutput= tree;  
178   AliTPCClustersRow clrow;
179   AliTPCClustersRow *pclrow=&clrow;  
180   clrow.SetClass("AliTPCclusterMI");
181   clrow.SetArray(1); // to make Clones array
182   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
183 }
184
185
186 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
187   // sigma y2 = in digits  - we don't know the angle
188   Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
189   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
190     (fPadWidth*fPadWidth);
191   Float_t sres = 0.25;
192   Float_t res = sd2+sres;
193   return res;
194 }
195
196
197 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
198   //sigma z2 = in digits - angle estimated supposing vertex constraint
199   Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
200   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
201   Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
202   angular*=angular;
203   angular/=12.;
204   Float_t sres = fParam->GetZSigma()/fZWidth;
205   sres *=sres;
206   Float_t res = angular +sd2+sres;
207   return res;
208 }
209
210 void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
211 AliTPCclusterMI &c) 
212 {
213   //
214   //  k    - Make cluster at position k  
215   //  bins - 2 D array of signals mapped to 1 dimensional array - 
216   //  max  - the number of time bins er one dimension
217   //  c    - refernce to cluster to be filled
218   //
219   Int_t i0=k/max;  //central pad
220   Int_t j0=k%max;  //central time bin
221
222   // set pointers to data
223   //Int_t dummy[5] ={0,0,0,0,0};
224   Float_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
225   for (Int_t di=-2;di<=2;di++){
226     matrix[di+2]  =  &bins[k+di*max];
227   }
228   //build matrix with virtual charge
229   Float_t sigmay2= GetSigmaY2(j0);
230   Float_t sigmaz2= GetSigmaZ2(j0);
231
232   Float_t vmatrix[5][5];
233   vmatrix[2][2] = matrix[2][0];
234   c.SetType(0);
235   c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
236   for (Int_t di =-1;di <=1;di++)
237     for (Int_t dj =-1;dj <=1;dj++){
238       Float_t amp = matrix[di+2][dj];
239       if ( (amp<2) && (fLoop<2)){
240         // if under threshold  - calculate virtual charge
241         Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
242         amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
243         if (amp>2) amp = 2;
244         vmatrix[2+di][2+dj]=amp;
245         vmatrix[2+2*di][2+2*dj]=0;
246         if ( (di*dj)!=0){       
247           //DIAGONAL ELEMENTS
248           vmatrix[2+2*di][2+dj] =0;
249           vmatrix[2+di][2+2*dj] =0;
250         }
251         continue;
252       }
253       if (amp<4){
254         //if small  amplitude - below  2 x threshold  - don't consider other one        
255         vmatrix[2+di][2+dj]=amp;
256         vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
257         if ( (di*dj)!=0){       
258           //DIAGONAL ELEMENTS
259           vmatrix[2+2*di][2+dj] =0;
260           vmatrix[2+di][2+2*dj] =0;
261         }
262         continue;
263       }
264       //if bigger then take everything
265       vmatrix[2+di][2+dj]=amp;
266       vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;      
267       if ( (di*dj)!=0){       
268           //DIAGONAL ELEMENTS
269           vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
270           vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
271         }      
272     }
273
274
275   
276   Float_t sumw=0;
277   Float_t sumiw=0;
278   Float_t sumi2w=0;
279   Float_t sumjw=0;
280   Float_t sumj2w=0;
281   //
282   for (Int_t i=-2;i<=2;i++)
283     for (Int_t j=-2;j<=2;j++){
284       Float_t amp = vmatrix[i+2][j+2];
285
286       sumw    += amp;
287       sumiw   += i*amp;
288       sumi2w  += i*i*amp;
289       sumjw   += j*amp;
290       sumj2w  += j*j*amp;
291     }    
292   //   
293   Float_t meani = sumiw/sumw;
294   Float_t mi2   = sumi2w/sumw-meani*meani;
295   Float_t meanj = sumjw/sumw;
296   Float_t mj2   = sumj2w/sumw-meanj*meanj;
297   //
298   Float_t ry = mi2/sigmay2;
299   Float_t rz = mj2/sigmaz2;
300   
301   //
302   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
303   if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) {
304     //
305     //if cluster looks like expected or Unfolding not switched on
306     //standard COG is used
307     //+1.2 deviation from expected sigma accepted
308     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
309
310     meani +=i0;
311     meanj +=j0;
312     //set cluster parameters
313     c.SetQ(sumw);
314     c.SetPad(meani-2.5);
315     c.SetTimeBin(meanj-3);
316     c.SetSigmaY2(mi2);
317     c.SetSigmaZ2(mj2);
318     c.SetType(0);
319     AddCluster(c,(Float_t*)vmatrix,k);
320     return;     
321   }
322   //
323   //unfolding when neccessary  
324   //
325   
326   Float_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
327   Float_t dummy[7]={0,0,0,0,0,0};
328   for (Int_t di=-3;di<=3;di++){
329     matrix2[di+3] =  &bins[k+di*max];
330     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
331     if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
332   }
333   Float_t vmatrix2[5][5];
334   Float_t sumu;
335   Float_t overlap;
336   UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
337   //
338   //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
339   meani +=i0;
340   meanj +=j0;
341   //set cluster parameters
342   c.SetQ(sumu);
343   c.SetPad(meani-2.5);
344   c.SetTimeBin(meanj-3);
345   c.SetSigmaY2(mi2);
346   c.SetSigmaZ2(mj2);
347   c.SetType(Char_t(overlap)+1);
348   AddCluster(c,(Float_t*)vmatrix,k);
349
350   //unfolding 2
351   meani-=i0;
352   meanj-=j0;
353   if (gDebug>4)
354     printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
355 }
356
357
358
359 void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, 
360                                       Float_t & sumu, Float_t & overlap )
361 {
362   //
363   //unfold cluster from input matrix
364   //data corresponding to cluster writen in recmatrix
365   //output meani and meanj
366
367   //take separatelly y and z
368
369   Float_t sum3i[7] = {0,0,0,0,0,0,0};
370   Float_t sum3j[7] = {0,0,0,0,0,0,0};
371
372   for (Int_t k =0;k<7;k++)
373     for (Int_t l = -1; l<=1;l++){
374       sum3i[k]+=matrix2[k][l];
375       sum3j[k]+=matrix2[l+3][k-3];
376     }
377   Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
378   //
379   //unfold  y 
380   Float_t sum3wi    = 0;  //charge minus overlap
381   Float_t sum3wio   = 0;  //full charge
382   Float_t sum3iw    = 0;  //sum for mean value
383   for (Int_t dk=-1;dk<=1;dk++){
384     sum3wio+=sum3i[dk+3];
385     if (dk==0){
386       sum3wi+=sum3i[dk+3];     
387     }
388     else{
389       Float_t ratio =1;
390       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
391            sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
392         Float_t xm2 = sum3i[-dk+3];
393         Float_t xm1 = sum3i[+3];
394         Float_t x1  = sum3i[2*dk+3];
395         Float_t x2  = sum3i[3*dk+3];    
396         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);      
397         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
398         ratio = w11/(w11+w12);   
399         for (Int_t dl=-1;dl<=1;dl++)
400           mratio[dk+1][dl+1] *= ratio;
401       }
402       Float_t amp = sum3i[dk+3]*ratio;
403       sum3wi+=amp;
404       sum3iw+= dk*amp;      
405     }
406   }
407   meani = sum3iw/sum3wi;
408   Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
409
410
411
412   //unfold  z 
413   Float_t sum3wj    = 0;  //charge minus overlap
414   Float_t sum3wjo   = 0;  //full charge
415   Float_t sum3jw    = 0;  //sum for mean value
416   for (Int_t dk=-1;dk<=1;dk++){
417     sum3wjo+=sum3j[dk+3];
418     if (dk==0){
419       sum3wj+=sum3j[dk+3];     
420     }
421     else{
422       Float_t ratio =1;
423       if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
424            (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
425         Float_t xm2 = sum3j[-dk+3];
426         Float_t xm1 = sum3j[+3];
427         Float_t x1  = sum3j[2*dk+3];
428         Float_t x2  = sum3j[3*dk+3];    
429         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);      
430         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
431         ratio = w11/(w11+w12);   
432         for (Int_t dl=-1;dl<=1;dl++)
433           mratio[dl+1][dk+1] *= ratio;
434       }
435       Float_t amp = sum3j[dk+3]*ratio;
436       sum3wj+=amp;
437       sum3jw+= dk*amp;      
438     }
439   }
440   meanj = sum3jw/sum3wj;
441   Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;  
442   overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);  
443   sumu = (sum3wj+sum3wi)/2.;
444   
445   if (overlap ==3) {
446     //if not overlap detected remove everything
447     for (Int_t di =-2; di<=2;di++)
448       for (Int_t dj =-2; dj<=2;dj++){
449         recmatrix[di+2][dj+2] = matrix2[3+di][dj];
450       }
451   }
452   else{
453     for (Int_t di =-1; di<=1;di++)
454       for (Int_t dj =-1; dj<=1;dj++){
455         Float_t ratio =1;
456         if (mratio[di+1][dj+1]==1){
457           recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
458           if (TMath::Abs(di)+TMath::Abs(dj)>1){
459             recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
460             recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
461           }       
462           recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
463         }
464         else
465           {
466             //if we have overlap in direction
467             recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];    
468             if (TMath::Abs(di)+TMath::Abs(dj)>1){
469               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
470               recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
471               //
472               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
473               recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
474             }
475             else{
476               ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
477               recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
478             }
479           }
480       }
481   }
482   if (gDebug>4) 
483     printf("%f\n", recmatrix[2][2]);
484   
485 }
486
487 Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
488 {
489   //
490   // estimate max
491   Float_t sumteor= 0;
492   Float_t sumamp = 0;
493
494   for (Int_t di = -1;di<=1;di++)
495     for (Int_t dj = -1;dj<=1;dj++){
496       if (vmatrix[2+di][2+dj]>2){
497         Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
498         sumteor += teor*vmatrix[2+di][2+dj];
499         sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
500       }
501     }    
502   Float_t max = sumamp/sumteor;
503   return max;
504 }
505
506 void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
507   //
508   //
509   // Transform cluster to the rotated global coordinata
510   // Assign labels to the cluster
511   // add the cluster to the array
512   // for more details - See  AliTPCTranform::Transform(x,i,0,1) 
513   Float_t meani = c.GetPad();
514   Float_t meanj = c.GetTimeBin();
515
516   Int_t ki = TMath::Nint(meani);
517   if (ki<0) ki=0;
518   if (ki>=fMaxPad) ki = fMaxPad-1;
519   Int_t kj = TMath::Nint(meanj);
520   if (kj<0) kj=0;
521   if (kj>=fMaxTime-3) kj=fMaxTime-4;
522   // ki and kj shifted as integers coordinata
523   if (fRowDig) {
524     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
525     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
526     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
527   }
528   
529   c.SetRow(fRow);
530   c.SetDetector(fSector);
531   Float_t s2 = c.GetSigmaY2();
532   Float_t w=fParam->GetPadPitchWidth(fSector);
533   c.SetSigmaY2(s2*w*w);
534   s2 = c.GetSigmaZ2(); 
535   c.SetSigmaZ2(s2*fZWidth*fZWidth);
536   //
537   //
538   //
539   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
540   if (!transform) {
541     AliFatal("Tranformations not in calibDB");
542   }
543   Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
544   Int_t i[1]={fSector};
545   transform->Transform(x,i,0,1);
546   c.SetX(x[0]);
547   c.SetY(x[1]);
548   c.SetZ(x[2]);
549   //
550   //
551   if (!fRecoParam->GetBYMirror()){
552     if (fSector%36>17){
553       c.SetY(-c.GetY());
554     }
555   }
556
557   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
558     c.SetType(-(c.GetType()+3));  //edge clusters
559   }
560   if (fLoop==2) c.SetType(100);
561
562   TClonesArray * arr = fRowCl->GetArray();
563   AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
564   if (fRecoParam->DumpSignal() &&matrix ) {
565     Int_t nbins=0;
566     Float_t *graph =0;
567     if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
568       nbins = fMaxTime;
569       graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
570     }
571     AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
572     cl->SetInfo(info);
573   }
574   if (!fRecoParam->DumpSignal()) {
575     cl->SetInfo(0);
576   }
577
578   fNcluster++;
579 }
580
581
582 //_____________________________________________________________________________
583 void AliTPCclustererMI::Digits2Clusters()
584 {
585   //-----------------------------------------------------------------
586   // This is a simple cluster finder.
587   //-----------------------------------------------------------------
588
589   if (!fInput) { 
590     Error("Digits2Clusters", "input tree not initialised");
591     return;
592   }
593  
594   if (!fOutput) {
595     Error("Digits2Clusters", "output tree not initialised");
596     return;
597   }
598
599   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
600   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
601   AliSimDigits digarr, *dummy=&digarr;
602   fRowDig = dummy;
603   fInput->GetBranch("Segment")->SetAddress(&dummy);
604   Stat_t nentries = fInput->GetEntries();
605   
606   fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3   after
607     
608   Int_t nclusters  = 0;
609
610   for (Int_t n=0; n<nentries; n++) {
611     fInput->GetEvent(n);
612     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
613       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
614       continue;
615     }
616     Int_t row = fRow;
617     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
618     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
619     //
620     AliTPCClustersRow *clrow= new AliTPCClustersRow();
621     fRowCl = clrow;
622     clrow->SetClass("AliTPCclusterMI");
623     clrow->SetArray(1);
624
625     clrow->SetID(digarr.GetID());
626     fOutput->GetBranch("Segment")->SetAddress(&clrow);
627     fRx=fParam->GetPadRowRadii(fSector,row);
628     
629     
630     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
631     fZWidth = fParam->GetZWidth();
632     if (fSector < kNIS) {
633       fMaxPad = fParam->GetNPadsLow(row);
634       fSign = (fSector < kNIS/2) ? 1 : -1;
635       fPadLength = fParam->GetPadPitchLength(fSector,row);
636       fPadWidth = fParam->GetPadPitchWidth();
637     } else {
638       fMaxPad = fParam->GetNPadsUp(row);
639       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
640       fPadLength = fParam->GetPadPitchLength(fSector,row);
641       fPadWidth  = fParam->GetPadPitchWidth();
642     }
643     
644     
645     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
646     fBins    =new Float_t[fMaxBin];
647     fSigBins =new Int_t[fMaxBin];
648     fNSigBins = 0;
649     memset(fBins,0,sizeof(Float_t)*fMaxBin);
650     
651     if (digarr.First()) //MI change
652       do {
653         Float_t dig=digarr.CurrentDigit();
654         if (dig<=fParam->GetZeroSup()) continue;
655         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
656         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
657         Int_t bin = i*fMaxTime+j;
658         fBins[bin]=dig/gain;
659         fSigBins[fNSigBins++]=bin;
660       } while (digarr.Next());
661     digarr.ExpandTrackBuffer();
662
663     FindClusters(noiseROC);
664
665     fOutput->Fill();
666     delete clrow;    
667     nclusters+=fNcluster;    
668     delete[] fBins;
669     delete[] fSigBins;
670   }  
671
672   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
673 }
674
675 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
676 {
677 //-----------------------------------------------------------------
678 // This is a cluster finder for the TPC raw data.
679 // The method assumes NO ordering of the altro channels.
680 // The pedestal subtraction can be switched on and off
681 // using an option of the TPC reconstructor
682 //-----------------------------------------------------------------
683
684   if (!fOutput) {
685     Error("Digits2Clusters", "output tree not initialised");
686     return;
687   }
688
689   fRowDig = NULL;
690   AliTPCROC * roc = AliTPCROC::Instance();
691   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
692   AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
693   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
694   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
695   //
696   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
697   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
698   if (fEventHeader){
699     fTimeStamp = fEventHeader->Get("Timestamp");  
700     fEventType = fEventHeader->Get("Type");  
701   }
702
703
704   Int_t nclusters  = 0;
705   
706   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
707   const Int_t kNIS = fParam->GetNInnerSector();
708   const Int_t kNOS = fParam->GetNOuterSector();
709   const Int_t kNS = kNIS + kNOS;
710   fZWidth = fParam->GetZWidth();
711   Int_t zeroSup = fParam->GetZeroSup();
712   //
713   //alocate memory for sector - maximal case
714   //
715   Float_t** allBins = NULL;
716   Int_t** allSigBins = NULL;
717   Int_t*  allNSigBins = NULL;
718   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
719   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
720   allBins = new Float_t*[nRowsMax];
721   allSigBins = new Int_t*[nRowsMax];
722   allNSigBins = new Int_t[nRowsMax];
723   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
724     //
725     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
726     allBins[iRow] = new Float_t[maxBin];
727     memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
728     allSigBins[iRow] = new Int_t[maxBin];
729     allNSigBins[iRow]=0;
730   }
731   //
732   // Loop over sectors
733   //
734   for(fSector = 0; fSector < kNS; fSector++) {
735
736     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
737     AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
738     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
739     //check the presence of the calibration
740     if (!noiseROC ||!pedestalROC ) {
741       AliError(Form("Missing calibration per sector\t%d\n",fSector));
742       continue;
743     }
744     Int_t nRows = 0;
745     Int_t nDDLs = 0, indexDDL = 0;
746     if (fSector < kNIS) {
747       nRows = fParam->GetNRowLow();
748       fSign = (fSector < kNIS/2) ? 1 : -1;
749       nDDLs = 2;
750       indexDDL = fSector * 2;
751     }
752     else {
753       nRows = fParam->GetNRowUp();
754       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
755       nDDLs = 4;
756       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
757     }
758
759     for (Int_t iRow = 0; iRow < nRows; iRow++) {
760       Int_t maxPad;
761       if (fSector < kNIS)
762         maxPad = fParam->GetNPadsLow(iRow);
763       else
764         maxPad = fParam->GetNPadsUp(iRow);
765       
766       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
767       memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
768       allNSigBins[iRow] = 0;
769     }
770     
771     // Loas the raw data for corresponding DDLs
772     rawReader->Reset();
773     input.SetOldRCUFormat(fIsOldRCUFormat);
774     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
775     Int_t digCounter=0;
776     // Begin loop over altro data
777     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
778     Float_t gain =1;
779     Int_t lastPad=-1;
780     while (input.Next()) {
781       if (input.GetSector() != fSector)
782         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
783
784       
785       Int_t iRow = input.GetRow();
786       if (iRow < 0 || iRow >= nRows){
787         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
788                       iRow, 0, nRows -1));
789         continue;
790       }
791       //pad
792       Int_t iPad = input.GetPad();
793       if (iPad < 0 || iPad >= nPadsMax) {
794         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
795                       iPad, 0, nPadsMax-1));
796         continue;
797       }
798       if (iPad!=lastPad){
799         gain    = gainROC->GetValue(iRow,iPad);
800         lastPad = iPad;
801       }
802       iPad+=3;
803       //time
804       Int_t iTimeBin = input.GetTime();
805       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
806         continue;
807         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
808                       iTimeBin, 0, iTimeBin -1));
809       }
810       iTimeBin+=3;
811
812       //signal
813       Float_t signal = input.GetSignal();
814       if (!calcPedestal && signal <= zeroSup) continue;      
815       if (!calcPedestal) {
816         Int_t bin = iPad*fMaxTime+iTimeBin;
817         allBins[iRow][bin] = signal/gain;
818         allSigBins[iRow][allNSigBins[iRow]++] = bin;
819       }else{
820         allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
821       }
822       allBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
823
824       // Temporary
825       digCounter++;
826     } // End of the loop over altro data
827     //
828     //
829     //
830     //
831     // Now loop over rows and perform pedestal subtraction
832     if (digCounter==0) continue;
833     //    if (calcPedestal) {
834     if (kTRUE) {
835       for (Int_t iRow = 0; iRow < nRows; iRow++) {
836         Int_t maxPad;
837         if (fSector < kNIS)
838           maxPad = fParam->GetNPadsLow(iRow);
839         else
840           maxPad = fParam->GetNPadsUp(iRow);
841
842         for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
843           //
844           // Temporary fix for data production - !!!! MARIAN
845           // The noise calibration should take mean and RMS - currently the Gaussian fit used
846           // In case of double peak  - the pad should be rejected
847           //
848           // Line mean - if more than given digits over threshold - make a noise calculation
849           // and pedestal substration
850           if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
851           //
852           if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
853           Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
854           //Float_t pedestal = TMath::Median(fMaxTime, p);      
855           Int_t id[3] = {fSector, iRow, iPad-3};
856           // calib values
857           Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
858           Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
859           Double_t rmsEvent       = rmsCalib;
860           Double_t pedestalEvent  = pedestalCalib;
861           ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
862           if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
863           if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;  
864           
865           //
866           for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
867             Int_t bin = iPad*fMaxTime+iTimeBin;
868             allBins[iRow][bin] -= pedestalEvent;
869             if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())  
870               allBins[iRow][bin] = 0;
871             if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())  
872               allBins[iRow][bin] = 0;
873             if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
874               allBins[iRow][bin] = 0;
875             if (allBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
876               allBins[iRow][bin] = 0;
877             if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
878           }
879         }
880       }
881     }
882     // Now loop over rows and find clusters
883     for (fRow = 0; fRow < nRows; fRow++) {
884       fRowCl = new AliTPCClustersRow;
885       fRowCl->SetClass("AliTPCclusterMI");
886       fRowCl->SetArray(1);
887       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
888       fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
889
890       fRx = fParam->GetPadRowRadii(fSector, fRow);
891       fPadLength = fParam->GetPadPitchLength(fSector, fRow);
892       fPadWidth  = fParam->GetPadPitchWidth();
893       if (fSector < kNIS)
894         fMaxPad = fParam->GetNPadsLow(fRow);
895       else
896         fMaxPad = fParam->GetNPadsUp(fRow);
897       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
898
899       fBins = allBins[fRow];
900       fSigBins = allSigBins[fRow];
901       fNSigBins = allNSigBins[fRow];
902
903       FindClusters(noiseROC);
904
905       fOutput->Fill();
906       delete fRowCl;    
907       nclusters += fNcluster;    
908     } // End of loop to find clusters
909   } // End of loop over sectors
910   
911   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
912     delete [] allBins[iRow];
913     delete [] allSigBins[iRow];
914   }  
915   delete [] allBins;
916   delete [] allSigBins;
917   delete [] allNSigBins;
918   
919   if (rawReader->GetEventId() && fOutput ){
920     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
921   }
922
923 }
924
925 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
926 {
927   
928   //
929   // add virtual charge at the edge   
930   //
931   Double_t kMaxDumpSize = 500000;
932   if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
933   //
934   fNcluster=0;
935   fLoop=1;
936   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
937   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
938   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
939   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
940   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
941   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
942   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
943   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
944     Int_t i = fSigBins[iSig];
945     if (i%fMaxTime<=crtime) continue;
946     Float_t *b = &fBins[i];
947     //absolute custs
948     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima  
949     //
950     if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters 
951     if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
952     if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
953     //
954     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF) 
955     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)    
956     if (!IsMaximum(*b,fMaxTime,b)) continue;
957     //
958     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
959     if (noise>fRecoParam->GetMaxNoise()) continue;
960     // sigma cuts
961     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima  
962     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF 
963     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)    
964   
965     AliTPCclusterMI c(kFALSE);   // default cosntruction  without info
966     Int_t dummy=0;
967     MakeCluster(i, fMaxTime, fBins, dummy,c);
968     
969     //}
970   }
971 }
972
973
974 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
975   //
976   // process signal on given pad - + streaming of additional information in special mode
977   //
978   // id[0] - sector
979   // id[1] - row
980   // id[2] - pad 
981
982   //
983   // ESTIMATE pedestal and the noise
984   // 
985   const Int_t kPedMax = 100;
986   Double_t kMaxDebugSize = 5000000.;
987   Float_t  max    =  0;
988   Float_t  maxPos =  0;
989   Int_t    median =  -1;
990   Int_t    count0 =  0;
991   Int_t    count1 =  0;
992   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
993   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
994   Int_t    firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
995   //
996   UShort_t histo[kPedMax];
997   memset(histo,0,kPedMax*sizeof(UShort_t));
998   for (Int_t i=0; i<fMaxTime; i++){
999     if (signal[i]<=0) continue;
1000     if (signal[i]>max && i>firstBin) {
1001       max = signal[i];
1002       maxPos = i;
1003     }
1004     if (signal[i]>kPedMax-1) continue;
1005     histo[int(signal[i]+0.5)]++;
1006     count0++;
1007   }
1008   //
1009   for (Int_t i=1; i<kPedMax; i++){
1010     if (count1<count0*0.5) median=i;
1011     count1+=histo[i];
1012   }
1013   // truncated mean  
1014   //
1015   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
1016   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
1017   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
1018   //
1019   for (Int_t idelta=1; idelta<10; idelta++){
1020     if (median-idelta<=0) continue;
1021     if (median+idelta>kPedMax) continue;
1022     if (count06<0.6*count1){
1023       count06+=histo[median-idelta];
1024       mean06 +=histo[median-idelta]*(median-idelta);
1025       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1026       count06+=histo[median+idelta];
1027       mean06 +=histo[median+idelta]*(median+idelta);
1028       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1029     }
1030     if (count09<0.9*count1){
1031       count09+=histo[median-idelta];
1032       mean09 +=histo[median-idelta]*(median-idelta);
1033       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1034       count09+=histo[median+idelta];
1035       mean09 +=histo[median+idelta]*(median+idelta);
1036       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1037     }
1038     if (count10<0.95*count1){
1039       count10+=histo[median-idelta];
1040       mean +=histo[median-idelta]*(median-idelta);
1041       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1042       count10+=histo[median+idelta];
1043       mean +=histo[median+idelta]*(median+idelta);
1044       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1045     }
1046   }
1047   if (count10) {
1048     mean  /=count10;
1049     rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1050   }
1051   if (count06) {
1052     mean06/=count06;
1053     rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1054   }
1055   if (count09) {
1056     mean09/=count09;
1057     rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1058   }
1059   rmsEvent = rms09;
1060   //
1061   pedestalEvent = median;
1062   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1063   //
1064   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1065   //
1066   // Dump mean signal info
1067   //
1068   (*fDebugStreamer)<<"Signal"<<
1069     "TimeStamp="<<fTimeStamp<<
1070     "EventType="<<fEventType<<
1071     "Sector="<<uid[0]<<
1072     "Row="<<uid[1]<<
1073     "Pad="<<uid[2]<<
1074     "Max="<<max<<
1075     "MaxPos="<<maxPos<<
1076     //
1077     "Median="<<median<<
1078     "Mean="<<mean<<
1079     "RMS="<<rms<<      
1080     "Mean06="<<mean06<<
1081     "RMS06="<<rms06<<
1082     "Mean09="<<mean09<<
1083     "RMS09="<<rms09<<
1084     "RMSCalib="<<rmsCalib<<
1085     "PedCalib="<<pedestalCalib<<
1086     "\n";
1087   //
1088   // fill pedestal histogram
1089   //
1090   AliTPCROC * roc = AliTPCROC::Instance();
1091
1092   //
1093   //
1094   //
1095   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1096   Float_t *dsignal = new Float_t[nchannels];
1097   Float_t *dtime   = new Float_t[nchannels];
1098   for (Int_t i=0; i<nchannels; i++){
1099     dtime[i] = i;
1100     dsignal[i] = signal[i];
1101   }
1102
1103   TGraph * graph;
1104   //
1105   // Big signals dumping
1106   //    
1107   if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin()) 
1108     (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1109       "TimeStamp="<<fTimeStamp<<
1110       "EventType="<<fEventType<<
1111       "Sector="<<uid[0]<<
1112       "Row="<<uid[1]<<
1113       "Pad="<<uid[2]<<
1114       "Graph="<<graph<<
1115       "Max="<<max<<
1116       "MaxPos="<<maxPos<<       
1117       //
1118       "Median="<<median<<
1119       "Mean="<<mean<<
1120       "RMS="<<rms<<      
1121       "Mean06="<<mean06<<
1122       "RMS06="<<rms06<<
1123       "Mean09="<<mean09<<
1124       "RMS09="<<rms09<<
1125       "\n";
1126   delete graph;  
1127
1128   delete [] dsignal;
1129   delete [] dtime;
1130   if (rms06>fRecoParam->GetMaxNoise()) {
1131     pedestalEvent+=1024.;
1132     return 1024+median; // sign noisy channel in debug mode
1133   }
1134   return median;
1135 }
1136
1137
1138
1139
1140