Update of tracking code provided by Sergei
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV0.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 $Log$
18 Revision 1.10  2001/12/05 15:04:34  hristov
19 Changes related to the corrections of AliRecPoint
20
21 Revision 1.9  2001/05/28 17:07:58  hristov
22 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
23
24 Revision 1.8  2001/05/07 08:06:44  cblume
25 Speedup of the code. Create only AliTRDcluster
26
27 Revision 1.7  2001/02/14 18:22:26  cblume
28 Change in the geometry of the padplane
29
30 Revision 1.6  2000/11/01 14:53:20  cblume
31 Merge with TRD-develop
32
33 Revision 1.1.4.6  2000/10/16 01:16:53  cblume
34 Changed timebin 0 to be the one closest to the readout
35
36 Revision 1.1.4.5  2000/10/15 23:40:01  cblume
37 Remove AliTRDconst
38
39 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
40 Made Getters const
41
42 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
43 Replace include files by forward declarations
44
45 Revision 1.1.4.2  2000/09/22 14:49:49  cblume
46 Adapted to tracking code
47
48 Revision 1.5  2000/10/02 21:28:19  fca
49 Removal of useless dependecies via forward declarations
50
51 Revision 1.4  2000/06/08 18:32:58  cblume
52 Make code compliant to coding conventions
53
54 Revision 1.3  2000/06/07 16:27:01  cblume
55 Try to remove compiler warnings on Sun and HP
56
57 Revision 1.2  2000/05/08 16:17:27  cblume
58 Merge TRD-develop
59
60 Revision 1.1.4.1  2000/05/08 15:08:41  cblume
61 Replace AliTRDcluster by AliTRDrecPoint
62
63 Revision 1.4  2000/06/08 18:32:58  cblume
64 Make code compliant to coding conventions
65
66 Revision 1.3  2000/06/07 16:27:01  cblume
67 Try to remove compiler warnings on Sun and HP
68
69 Revision 1.2  2000/05/08 16:17:27  cblume
70 Merge TRD-develop
71
72 Revision 1.1.4.1  2000/05/08 15:08:41  cblume
73 Replace AliTRDcluster by AliTRDrecPoint
74
75 Revision 1.1  2000/02/28 18:58:33  cblume
76 Add new TRD classes
77
78 */
79
80 ///////////////////////////////////////////////////////////////////////////////
81 //                                                                           //
82 // TRD cluster finder for the fast simulator. It takes the hits from the     //
83 // fast simulator (one hit per plane) and transforms them                    //
84 // into cluster, by applying position smearing and merging                   //
85 // of nearby cluster. The searing is done uniformly in z-direction           //
86 // over the length of a readout pad. In rphi-direction a Gaussian            //
87 // smearing is applied with a sigma given by fRphiSigma.                     //
88 // Clusters are considered as overlapping when they are closer in            //
89 // rphi-direction than the value defined in fRphiDist.                       //
90 // Use the macro fastClusterCreate.C to create the cluster.                  //
91 //                                                                           //
92 ///////////////////////////////////////////////////////////////////////////////
93
94 #include <TRandom.h>
95 #include <TTree.h>
96  
97 #include "AliRun.h"
98
99 #include "AliTRD.h"
100 #include "AliTRDclusterizerV0.h"
101 #include "AliTRDhit.h"
102 #include "AliTRDgeometry.h"
103 #include "AliTRDrecPoint.h"
104 #include "AliTRDparameter.h"
105
106 ClassImp(AliTRDclusterizerV0)
107
108 //_____________________________________________________________________________
109 AliTRDclusterizerV0::AliTRDclusterizerV0():AliTRDclusterizer()
110 {
111   //
112   // AliTRDclusterizerV0 default constructor
113   //
114
115 }
116
117 //_____________________________________________________________________________
118 AliTRDclusterizerV0::AliTRDclusterizerV0(const Text_t* name, const Text_t* title)
119                     :AliTRDclusterizer(name,title)
120 {
121   //
122   // AliTRDclusterizerV0 default constructor
123   //
124
125   Init();
126
127 }
128
129 //_____________________________________________________________________________
130 AliTRDclusterizerV0::~AliTRDclusterizerV0()
131 {
132   //
133   // AliTRDclusterizerV0 destructor
134   //
135
136 }
137
138 //_____________________________________________________________________________
139 void AliTRDclusterizerV0::Init()
140 {
141   //
142   // Initializes the cluster finder
143   //
144
145   // Position resolution in rphi-direction
146   fRphiSigma  = 0.02;
147   // Minimum distance of non-overlapping cluster
148   fRphiDist   = 1.0;
149
150 }
151
152 //_____________________________________________________________________________
153 Bool_t AliTRDclusterizerV0::MakeClusters()
154 {
155   //
156   // Generates the cluster
157   //
158
159   if (fTRD->IsVersion() != 0) {
160     printf("AliTRDclusterizerV0::MakeCluster -- ");
161     printf("TRD must be version 0 (fast simulator).\n");
162     return kFALSE; 
163   }
164
165   // Get the geometry
166   AliTRDgeometry *geo = fTRD->GetGeometry();
167
168   // Create a default parameter class if none is defined
169   if (!fPar) {
170     fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
171     if (fVerbose > 0) {
172       printf("<AliTRDclusterizerV0::MakeCluster> ");
173       printf("Create the default parameter object.\n");
174     }
175   }
176   
177   printf("AliTRDclusterizerV0::MakeCluster -- ");
178   printf("Start creating cluster.\n");
179
180   Int_t nBytes = 0;
181
182   AliTRDhit *hit;
183   
184   // Get the pointer to the hit tree
185   TTree     *hitTree      = gAlice->TreeH();
186   // Get the pointer to the reconstruction tree
187   TTree     *clusterTree  = gAlice->TreeR();
188
189   TObjArray *chamberArray = new TObjArray();
190
191   // Get the number of entries in the hit tree
192   // (Number of primary particles creating a hit somewhere)
193   Int_t nTrack = (Int_t) hitTree->GetEntries();
194
195   // Loop through all the chambers
196   for (Int_t icham = 0; icham < AliTRDgeometry::Ncham(); icham++) {
197     for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) {
198       for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
199
200         Int_t   nColMax     = fPar->GetColMax(iplan);
201         Float_t row0        = fPar->GetRow0(iplan,icham,isect);
202         Float_t col0        = fPar->GetCol0(iplan);
203         Float_t time0       = fPar->GetTime0(iplan);
204
205         Float_t rowPadSize  = fPar->GetRowPadSize(iplan,icham,isect);
206         Float_t colPadSize  = fPar->GetColPadSize(iplan);
207         Float_t timeBinSize = fPar->GetTimeBinSize();
208
209         // Loop through all entries in the tree
210         for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
211
212           gAlice->ResetHits();
213           nBytes += hitTree->GetEvent(iTrack);
214
215           // Get the number of hits in the TRD created by this particle
216           Int_t nHit = fTRD->Hits()->GetEntriesFast();
217
218           // Loop through the TRD hits  
219           for (Int_t iHit = 0; iHit < nHit; iHit++) {
220
221             if (!(hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit))) 
222               continue;
223
224             Float_t pos[3];
225                     pos[0]   = hit->X();
226                     pos[1]   = hit->Y();
227                     pos[2]   = hit->Z();
228             Int_t   track    = hit->Track();
229             Int_t   detector = hit->GetDetector();
230             Int_t   plane    = geo->GetPlane(detector);
231             Int_t   sector   = geo->GetSector(detector);
232             Int_t   chamber  = geo->GetChamber(detector);        
233
234             if ((sector  != isect) ||
235                 (plane   != iplan) ||
236                 (chamber != icham)) 
237               continue;
238
239             // Rotate the sectors on top of each other
240             Float_t rot[3];
241             geo->Rotate(detector,pos,rot);
242
243             // Add this recPoint to the temporary array for this chamber
244             AliTRDrecPoint *recPoint = new AliTRDrecPoint("");
245             recPoint->SetLocalRow(rot[2]);
246             recPoint->SetLocalCol(rot[1]);
247             recPoint->SetLocalTime(rot[0]);
248             recPoint->SetEnergy(0);
249             recPoint->SetDetector(detector);
250             recPoint->AddDigit(track);
251             chamberArray->Add(recPoint);
252
253           }
254
255         }
256   
257         // Loop through the temporary cluster-array
258         for (Int_t iClus1 = 0; iClus1 < chamberArray->GetEntries(); iClus1++) {
259
260           AliTRDrecPoint *recPoint1 = (AliTRDrecPoint *) 
261                                       chamberArray->UncheckedAt(iClus1);
262           Float_t row1  = recPoint1->GetLocalRow();
263           Float_t col1  = recPoint1->GetLocalCol();
264           Float_t time1 = recPoint1->GetLocalTime();
265
266           if (recPoint1->GetEnergy() < 0) continue;        // Skip marked cluster  
267
268           const Int_t kNsave  = 5;
269           Int_t idxSave[kNsave];
270           Int_t iSave = 0;
271
272           const Int_t kNsaveTrack = 3;
273           Int_t tracks[kNsaveTrack];
274           tracks[0] = recPoint1->GetDigit(0);
275
276           // Check the other cluster to see, whether there are close ones
277           for (Int_t iClus2 = iClus1 + 1; iClus2 < chamberArray->GetEntries(); iClus2++) {
278
279             AliTRDrecPoint *recPoint2 = (AliTRDrecPoint *) 
280                                         chamberArray->UncheckedAt(iClus2);
281             Float_t row2 = recPoint2->GetLocalRow();
282             Float_t col2 = recPoint2->GetLocalCol();
283
284             if ((TMath::Abs(row1 - row2) < rowPadSize) ||
285                 (TMath::Abs(col1 - col2) <  fRphiDist)) {
286               if (iSave == kNsave) {
287                 printf("AliTRDclusterizerV0::MakeCluster -- ");
288                 printf("Boundary error: iSave = %d, kNsave = %d.\n"
289                       ,iSave,kNsave);
290               }
291               else {                              
292                 idxSave[iSave]  = iClus2;
293                 iSave++;
294                 if (iSave < kNsaveTrack) tracks[iSave] = recPoint2->GetDigit(0);
295               }
296             }
297           }
298      
299           // Merge close cluster
300           Float_t rowMerge = row1;
301           Float_t colMerge = col1;
302           if (iSave) {
303             for (Int_t iMerge = 0; iMerge < iSave; iMerge++) {
304               AliTRDrecPoint *recPoint2 =
305                 (AliTRDrecPoint *) chamberArray->UncheckedAt(idxSave[iMerge]);
306               rowMerge += recPoint2->GetLocalRow();
307               colMerge += recPoint2->GetLocalCol();
308               recPoint2->SetEnergy(-1);     // Mark merged cluster
309             }
310             rowMerge /= (iSave + 1);
311             colMerge /= (iSave + 1);
312           }
313
314           Float_t smear[3];
315
316           // The position smearing in row-direction (uniform over pad width)            
317           Int_t row = (Int_t) ((rowMerge - row0) / rowPadSize);
318           smear[0]  = (row + gRandom->Rndm()) * rowPadSize + row0;
319
320           // The position smearing in rphi-direction (Gaussian)
321           smear[1] = 0;
322           do
323             smear[1] = gRandom->Gaus(colMerge,fRphiSigma);
324           while ((smear[1] < col0                        ) ||
325                  (smear[1] > col0 + nColMax * colPadSize));
326
327           // Time direction stays unchanged
328           smear[2] = time1;
329          
330           // Transform into local coordinates
331           smear[0] = (Int_t) ((smear[0] -  row0) /  rowPadSize);
332           smear[1] = (Int_t) ((smear[1] -  col0) /  colPadSize);
333           smear[2] = (Int_t) ((time0 - smear[2]) / timeBinSize);
334
335           // Add the smeared cluster to the output array 
336           Int_t   detector  = recPoint1->GetDetector();
337           Int_t   tr[9]     = { -1   };
338           Float_t pos[3];
339           Float_t sigma[2]  = {  0.0 };
340           pos[0] = smear[1];
341           pos[1] = smear[0];
342           pos[2] = (time0 - smear[2]) / timeBinSize;
343           fTRD->AddCluster(pos,detector,0.0,tr,sigma,0);
344
345         }
346
347         // Clear the temporary cluster-array and delete the cluster
348         chamberArray->Delete();
349
350       }
351     }
352   }
353
354   printf("AliTRDclusterizerV0::MakeCluster -- ");
355   printf("Found %d points.\n",fTRD->RecPoints()->GetEntries());
356   printf("AliTRDclusterizerV0::MakeCluster -- ");
357   printf("Fill the cluster tree.\n");
358   clusterTree->Fill();
359
360   return kTRUE;
361
362 }