Latest Fixes from M.Bondila
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSDD.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
17 #include "AliITSClusterFinderSDD.h"
18 #include "AliITSMapA2.h"
19 #include "AliITSMapA1.h"
20 #include "AliITS.h"
21 #include "AliRun.h"
22
23
24
25 ClassImp(AliITSClusterFinderSDD)
26
27 //----------------------------------------------------------
28 AliITSClusterFinderSDD::AliITSClusterFinderSDD
29 (AliITSsegmentation *seg, AliITSresponse *response, TClonesArray *digits, TClonesArray *recp)   
30 {
31   // constructor
32     fSegmentation=seg;
33     fResponse=response;
34     fDigits=digits;
35     fClusters=recp;
36     fNclusters= fClusters->GetEntriesFast();
37     SetCutAmplitude();
38     SetDAnode();
39     SetDTime();
40     SetMinPeak();
41     SetNCells();
42     SetMap();
43 }
44
45 //_____________________________________________________________________________
46 AliITSClusterFinderSDD::AliITSClusterFinderSDD()
47 {
48   // constructor
49     fSegmentation=0;
50     fResponse=0;
51     fDigits=0;
52     fClusters=0;
53     fNclusters=0;
54     fMap=0;
55     fMapA2=0;
56     SetCutAmplitude();
57     SetDAnode();
58     SetDTime();
59     SetMinPeak();
60     SetNCells();
61
62 }
63
64 //_____________________________________________________________________________
65 AliITSClusterFinderSDD::~AliITSClusterFinderSDD()
66 {
67     // destructor
68
69     if(fMap) delete fMap;
70     if(fMapA2) delete fMapA2;
71
72 }
73 //__________________________________________________________________________
74 AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){
75   //     Copy Constructor 
76   if(&source == this) return;
77   this->fClusters = source.fClusters ;
78   this->fNclusters = source.fNclusters ;
79   this->fMap = source.fMap ;
80   this->fMapA2 = source.fMapA2 ;
81   this->fCutAmplitude = source.fCutAmplitude ;
82   this->fDAnode = source.fDAnode ;
83   this->fDTime = source.fDTime ;
84   this->fMinPeak = source.fMinPeak ;
85   this->fMinNCells = source.fMinNCells ;
86   return;
87 }
88
89 //_________________________________________________________________________
90 AliITSClusterFinderSDD& 
91   AliITSClusterFinderSDD::operator=(const AliITSClusterFinderSDD &source) {
92   //    Assignment operator
93   if(&source == this) return *this;
94   this->fClusters = source.fClusters ;
95   this->fNclusters = source.fNclusters ;
96   this->fMap = source.fMap ;
97   this->fMapA2 = source.fMapA2 ;
98   this->fCutAmplitude = source.fCutAmplitude ;
99   this->fDAnode = source.fDAnode ;
100   this->fDTime = source.fDTime ;
101   this->fMinPeak = source.fMinPeak ;
102   this->fMinNCells = source.fMinNCells ;
103   return *this;
104 }
105
106 //_____________________________________________________________________________
107 void AliITSClusterFinderSDD::SetMap()
108 {
109   // set map
110     if(!fMapA2) fMapA2=new AliITSMapA2(fSegmentation,fDigits,(double)fCutAmplitude);
111     if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
112
113 }
114 //_____________________________________________________________________________
115
116 void AliITSClusterFinderSDD::Find1DClusters()
117 {
118   // find 1D clusters
119
120     AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
121
122    // retrieve the parameters 
123    Int_t fNofMaps = fSegmentation->Npz();
124    Int_t fMaxNofSamples = fSegmentation->Npx();
125    Int_t fNofAnodes = fNofMaps/2;
126    Int_t dummy=0;
127    Float_t fTimeStep = fSegmentation->Dpx(dummy);
128    Float_t fSddLength = fSegmentation->Dx();
129    Float_t fDriftSpeed = fResponse->DriftSpeed();
130    
131    Float_t anodePitch = fSegmentation->Dpz(dummy);
132    // map the signal
133    fMapA2->FillMap();
134
135    // Piergiorgio's code - do not subtract baseline since we start
136    // from digits and do not duplicate arrays, i.e. use fMap instead
137    // of Float_t fadc[2*fNofAnodes][fMaxNofSamples];
138  
139
140   Int_t nofFoundClusters = 0;
141   Int_t i;
142   Float_t **dfadc = new Float_t*[fNofMaps];
143   for(i=0;i<fNofMaps;i++) dfadc[i] = new Float_t[fMaxNofSamples];
144   Float_t fadc, fadc1, fadc2;
145   Int_t j,k,idx,l,m;
146   for(j=0;j<2;j++) {
147     for(k=0;k<fNofAnodes;k++) {
148       idx = j*fNofAnodes+k;
149       // signal (fadc) & derivative (dfadc)
150       for(l=0; l<fMaxNofSamples; l++) {
151         fadc2=(Float_t)fMapA2->GetSignal(idx,l);
152         fadc1=(Float_t)fMapA2->GetSignal(idx,l-1);
153         if(l>0) dfadc[k][l-1] = fadc2-fadc1;
154       } // samples
155     } // anodes
156     
157     for(k=0;k<fNofAnodes;k++) {
158       idx = j*fNofAnodes+k;
159  
160       Int_t imax = 0;
161       Int_t imaxd = 0;
162       Int_t it=0;
163       while(it <= fMaxNofSamples-3) {
164         //      cout << "sample: " << it << endl;
165         
166         imax = it;
167         imaxd = it;
168         // maximum of signal      
169         
170         Float_t fadcmax = 0.;
171         Float_t dfadcmax = 0.;
172         Int_t lthrmina = 1;
173         //      if(it >= 60) lthrmina = 2;
174         //      if(it >= 100) lthrmina = 3;
175         Int_t lthrmint = 2;
176         //if(it >= 60) lthrmint = 3;
177         //if(it >= 100) lthrmint = 4;
178
179         Int_t lthra = 0;
180         Int_t lthrt = 0;
181         for(m=0;m<10;m++) {
182           Int_t id = it+m;
183           if(id>=fMaxNofSamples) break;
184           fadc=fMapA2->GetSignal(idx,id);
185           if(fadc > fadcmax) {
186             fadcmax = fadc;
187             if(fadc > 0) { lthra++; lthrt++; }
188             imax = id;
189           }
190           if(dfadc[k][id] > dfadcmax) {
191             dfadcmax = dfadc[k][id];
192             imaxd = id;
193           }
194         }
195         it = imaxd;
196         // skip if no signal over threshold
197         if(fMapA2->TestHit(idx,imax) == kEmpty) {it++; continue;}
198
199         if(k>0) {
200           if(fMapA2->TestHit(idx-1,imax) != kEmpty) lthra++;
201         }
202         if(k<fNofAnodes-1)
203           if(fMapA2->TestHit(idx+1,imax) != kEmpty) lthra++;
204
205         if(imax>0) {
206           if(fMapA2->TestHit(idx,imax-1) != kEmpty) lthrt++;
207         }
208         if(imax<fMaxNofSamples)
209           if(fMapA2->TestHit(idx,imax+1) != kEmpty) lthrt++;
210         
211         // cluster charge
212         Int_t tstart = it-1;
213         
214         Bool_t ilcl = 0;
215         if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
216         
217         if(ilcl) {
218           nofFoundClusters++;
219           Int_t tstop = tstart;
220           Float_t dfadcmin = 10000.;
221           Int_t ij;
222           for(ij=0; ij<10; ij++) {
223             if(dfadc[k][it+ij] < dfadcmin) {
224               tstop = it+ij+1;
225               dfadcmin = dfadc[k][it+ij];
226             }
227           }
228           
229           Float_t clusterCharge = 0.;
230           Float_t clusterAnode = k+0.5;
231           Float_t clusterTime = 0.;
232           Float_t clusterMult = 0.;
233           Float_t clusterPeakAmplitude = 0.;
234           Int_t its;
235           Float_t n, baseline;
236           fResponse->GetNoiseParam(n,baseline);
237           for(its=tstart; its<=tstop; its++) {
238             fadc=fMapA2->GetSignal(idx,its)-baseline;
239             clusterCharge += fadc;
240             if(fadc > clusterPeakAmplitude) clusterPeakAmplitude = fadc;
241             clusterTime += fadc*its;
242             clusterMult++;
243             if(its == tstop) {
244                clusterTime /= (clusterCharge/fTimeStep);   // ns
245                clusterCharge *= (fTimeStep/160.);          // keV
246                if(clusterTime > 58.2) clusterTime -= 58.2;   // ns
247               /*
248               else {
249                 cout << "Warning: cluster with negative time " << clusterTime << ", peak ampl.: " << clusterPeakAmplitude << ", mult: " << clusterMult << ", charge: " << clusterCharge << endl;
250                 cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
251               }
252               */
253               }
254           }
255           // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
256
257           Float_t clusteranodePath = (0.06 + clusterAnode - fNofAnodes/2)*anodePitch;
258           Float_t clusterDriftPath = clusterTime*fDriftSpeed;
259           clusterDriftPath = fSddLength-clusterDriftPath;
260
261           if(clusterCharge < 0.) break;
262
263           //printf("wing clusterMult clusterAnode clusterTime %d %f %f %f \n",j+1,clusterMult, clusterAnode, clusterTime);
264
265           AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
266           //fClusters->Add(point);
267           iTS->AddCluster(1,clust);
268           it = tstop;
269         } // ilcl
270         
271         it++;
272         
273       } // while (samples)
274     } // anodes
275   } // detectors (2)
276
277
278   fMapA2->ClearMap();
279   
280   for(i=0;i<fNofMaps;i++) delete[] dfadc[i];
281   delete [] dfadc;
282
283   return;
284
285 }
286
287 //_____________________________________________________________________________
288 void  AliITSClusterFinderSDD::GroupClusters()
289 {
290   // group clusters
291   Int_t dummy=0;
292   Float_t fTimeStep = fSegmentation->Dpx(dummy);
293
294
295   // get number of clusters for this module
296   Int_t nofClusters = fClusters->GetEntriesFast();
297   nofClusters -= fNclusters;
298
299   AliITSRawClusterSDD *clusterI;
300   AliITSRawClusterSDD *clusterJ;
301
302   Int_t *label = new Int_t [nofClusters];
303   Int_t i,j;
304   for(i=0; i<nofClusters; i++) label[i] = 0;
305   for(i=0; i<nofClusters; i++) { 
306     if(label[i] != 0) continue;
307     for(j=i+1; j<nofClusters; j++) { 
308       if(label[j] != 0) continue;
309       clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
310       clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
311       // 1.3 good
312       if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
313       if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
314       Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
315       if(!pair) continue;
316       //      clusterI->PrintInfo();
317       //      clusterJ->PrintInfo();
318       clusterI->Add(clusterJ);
319       label[j] = 1;
320       fClusters->RemoveAt(j);
321     } // J clusters  
322     label[i] = 1;
323   } // I clusters
324   fClusters->Compress();
325   
326   delete [] label;
327   return;
328
329 }
330
331 //_____________________________________________________________________________
332
333 void AliITSClusterFinderSDD::SelectClusters()
334 {
335   // get number of clusters for this module
336   Int_t nofClusters = fClusters->GetEntriesFast();
337   nofClusters -= fNclusters;
338
339   Int_t i;
340   for(i=0; i<nofClusters; i++) { 
341     AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
342     Int_t rmflg = 0;
343     Float_t wy = 0.;
344     if(clusterI->Anodes() != 0.) {
345       wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
346     }
347     Float_t amp = clusterI->PeakAmpl();
348     if(amp < fMinPeak) rmflg = 1;  
349     if(wy < fMinNCells) rmflg = 1;
350     if(rmflg) fClusters->RemoveAt(i);
351   } // I clusters
352   fClusters->Compress();
353   return;
354
355 }
356
357 //_____________________________________________________________________________
358
359 void AliITSClusterFinderSDD::GetRecPoints()
360 {
361   // get rec points
362
363   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
364
365   // get number of clusters for this module
366   Int_t nofClusters = fClusters->GetEntriesFast();
367   nofClusters -= fNclusters;
368
369   const Float_t kconvGeV = 1.e-6; // GeV -> KeV
370   const Float_t kconv = 1.0e-4; 
371   const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
372   const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
373
374   fMap->FillMap();
375
376   Int_t i, ix, iz;
377   for(i=0; i<nofClusters; i++) { 
378     AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
379     fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
380     AliITSdigitSDD *dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1);
381     AliITSRecPoint rnew;
382     rnew.SetX(clusterI->X());
383     rnew.SetZ(clusterI->Z());
384     rnew.SetQ(clusterI->Q());   // in KeV - should be ADC
385     rnew.SetdEdX(kconvGeV*clusterI->Q());
386     rnew.SetSigmaX2(kRMSx*kRMSx);
387     rnew.SetSigmaZ2(kRMSz*kRMSz);
388     rnew.fTracks[0]=dig->fTracks[0];
389     rnew.fTracks[1]=dig->fTracks[1];
390     rnew.fTracks[2]=dig->fTracks[2];
391     iTS->AddRecPoint(rnew);
392   } // I clusters
393
394   fMap->ClearMap();
395 }
396
397 //_____________________________________________________________________________
398
399 void AliITSClusterFinderSDD::FindRawClusters()
400 {
401   // find raw clusters
402     Find1DClusters();
403     GroupClusters();
404     SelectClusters();
405     GetRecPoints();
406 }
407
408
409
410