3b14832b6b1a925cdc3f3848ebfe700a9a4f52b1
[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.1.4.1  2000/05/08 15:08:41  cblume
19 Replace AliTRDcluster by AliTRDrecPoint
20
21 Revision 1.1  2000/02/28 18:58:33  cblume
22 Add new TRD classes
23
24 */
25
26 ///////////////////////////////////////////////////////////////////////////////
27 //                                                                           //
28 // TRD cluster finder for the fast simulator. It takes the hits from the     //
29 // fast simulator (one hit per plane) and transforms them                    //
30 // into cluster, by applying position smearing and merging                   //
31 // of nearby cluster. The searing is done uniformly in z-direction           //
32 // over the length of a readout pad. In rphi-direction a Gaussian            //
33 // smearing is applied with a sigma given by fRphiSigma.                     //
34 // Clusters are considered as overlapping when they are closer in            //
35 // rphi-direction than the value defined in fRphiDist.                       //
36 // Use the macro fastClusterCreate.C to create the cluster.                  //
37 //                                                                           //
38 ///////////////////////////////////////////////////////////////////////////////
39
40 #include <TRandom.h>
41
42 #include "AliTRDclusterizerV0.h"
43 #include "AliTRDconst.h"
44 #include "AliTRDgeometry.h"
45 #include "AliTRDrecPoint.h"
46
47 ClassImp(AliTRDclusterizerV0)
48
49 //_____________________________________________________________________________
50 AliTRDclusterizerV0::AliTRDclusterizerV0():AliTRDclusterizer()
51 {
52   //
53   // AliTRDclusterizerV0 default constructor
54   //
55
56 }
57
58 //_____________________________________________________________________________
59 AliTRDclusterizerV0::AliTRDclusterizerV0(const Text_t* name, const Text_t* title)
60                     :AliTRDclusterizer(name,title)
61 {
62   //
63   // AliTRDclusterizerV0 default constructor
64   //
65
66   Init();
67
68 }
69
70 //_____________________________________________________________________________
71 AliTRDclusterizerV0::~AliTRDclusterizerV0()
72 {
73
74 }
75
76 //_____________________________________________________________________________
77 void AliTRDclusterizerV0::Init()
78 {
79   //
80   // Initializes the cluster finder
81   //
82
83   // Position resolution in rphi-direction
84   fRphiSigma  = 0.02;
85   // Minimum distance of non-overlapping cluster
86   fRphiDist   = 1.0;
87
88 }
89
90 //_____________________________________________________________________________
91 Bool_t AliTRDclusterizerV0::MakeCluster()
92 {
93   //
94   // Generates the cluster
95   //
96
97   // Get the pointer to the detector class and check for version 1
98   AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
99   if (TRD->IsVersion() != 0) {
100     printf("AliTRDclusterizerV0::MakeCluster -- ");
101     printf("TRD must be version 0 (fast simulator).\n");
102     return kFALSE; 
103   }
104
105   // Get the geometry
106   AliTRDgeometry *Geo = TRD->GetGeometry();
107   
108   printf("AliTRDclusterizerV0::MakeCluster -- ");
109   printf("Start creating cluster.\n");
110
111   Int_t nBytes = 0;
112
113   AliTRDhit      *Hit;
114   
115   // Get the pointer to the hit tree
116   TTree *HitTree     = gAlice->TreeH();
117   // Get the pointer to the reconstruction tree
118   TTree *ClusterTree = gAlice->TreeR();
119
120   TObjArray *Chamber = new TObjArray();
121
122   // Get the number of entries in the hit tree
123   // (Number of primary particles creating a hit somewhere)
124   Int_t nTrack = (Int_t) HitTree->GetEntries();
125
126   // Loop through all the chambers
127   for (Int_t icham = 0; icham < kNcham; icham++) {
128     for (Int_t iplan = 0; iplan < kNplan; iplan++) {
129       for (Int_t isect = 0; isect < kNsect; isect++) {
130
131         Int_t   nColMax     = Geo->GetColMax(iplan);
132         Float_t row0        = Geo->GetRow0(iplan,icham,isect);
133         Float_t col0        = Geo->GetCol0(iplan);
134         Float_t time0       = Geo->GetTime0(iplan);
135
136         Float_t rowPadSize  = Geo->GetRowPadSize();
137         Float_t colPadSize  = Geo->GetColPadSize();
138         Float_t timeBinSize = Geo->GetTimeBinSize();
139
140         // Loop through all entries in the tree
141         for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
142
143           gAlice->ResetHits();
144           nBytes += HitTree->GetEvent(iTrack);
145
146           // Get the number of hits in the TRD created by this particle
147           Int_t nHit = TRD->Hits()->GetEntriesFast();
148
149           // Loop through the TRD hits  
150           for (Int_t iHit = 0; iHit < nHit; iHit++) {
151
152             if (!(Hit = (AliTRDhit *) TRD->Hits()->UncheckedAt(iHit))) 
153               continue;
154
155             Float_t pos[3];
156                     pos[0]   = Hit->fX;
157                     pos[1]   = Hit->fY;
158                     pos[2]   = Hit->fZ;
159             Int_t   track    = Hit->fTrack;
160             Int_t   detector = Hit->fDetector;
161             Int_t   plane    = Geo->GetPlane(detector);
162             Int_t   sector   = Geo->GetSector(detector);
163             Int_t   chamber  = Geo->GetChamber(detector);        
164
165             if ((sector  != isect) ||
166                 (plane   != iplan) ||
167                 (chamber != icham)) 
168               continue;
169
170             // Rotate the sectors on top of each other
171             Float_t rot[3];
172             Geo->Rotate(detector,pos,rot);
173
174             // Add this recPoint to the temporary array for this chamber
175             AliTRDrecPoint *RecPoint = new AliTRDrecPoint();
176             RecPoint->SetLocalRow(rot[2]);
177             RecPoint->SetLocalCol(rot[1]);
178             RecPoint->SetLocalTime(rot[0]);
179             RecPoint->SetEnergy(0);
180             RecPoint->SetDetector(detector);
181             RecPoint->AddDigit(track);
182             Chamber->Add(RecPoint);
183
184           }
185
186         }
187   
188         // Loop through the temporary cluster-array
189         for (Int_t iClus1 = 0; iClus1 < Chamber->GetEntries(); iClus1++) {
190
191           AliTRDrecPoint *RecPoint1 = (AliTRDrecPoint *) Chamber->UncheckedAt(iClus1);
192           Float_t row1  = RecPoint1->GetLocalRow();
193           Float_t col1  = RecPoint1->GetLocalCol();
194           Float_t time1 = RecPoint1->GetLocalTime();
195
196           if (RecPoint1->GetEnergy() < 0) continue;        // Skip marked cluster  
197
198           const Int_t nSave  = 5;
199           Int_t idxSave[nSave];
200           Int_t iSave = 0;
201
202           const Int_t nTrack = 3;
203           Int_t tracks[nTrack];
204           tracks[0] = RecPoint1->GetDigit(0);
205
206           // Check the other cluster to see, whether there are close ones
207           for (Int_t iClus2 = iClus1 + 1; iClus2 < Chamber->GetEntries(); iClus2++) {
208
209             AliTRDrecPoint *RecPoint2 = (AliTRDrecPoint *) Chamber->UncheckedAt(iClus2);
210             Float_t row2 = RecPoint2->GetLocalRow();
211             Float_t col2 = RecPoint2->GetLocalCol();
212
213             if ((TMath::Abs(row1 - row2) < rowPadSize) ||
214                 (TMath::Abs(col1 - col2) <  fRphiDist)) {
215               if (iSave == nSave) {
216                 printf("AliTRDclusterizerV0::MakeCluster -- ");
217                 printf("Boundary error: iSave = %d, nSave = %d.\n"
218                       ,iSave,nSave);
219               }
220               else {                              
221                 idxSave[iSave]  = iClus2;
222                 iSave++;
223                 if (iSave < nTrack) tracks[iSave] = RecPoint2->GetDigit(0);
224               }
225             }
226           }
227      
228           // Merge close cluster
229           Float_t rowMerge = row1;
230           Float_t colMerge = col1;
231           if (iSave) {
232             for (Int_t iMerge = 0; iMerge < iSave; iMerge++) {
233               AliTRDrecPoint *RecPoint2 =
234                 (AliTRDrecPoint *) Chamber->UncheckedAt(idxSave[iMerge]);
235               rowMerge += RecPoint2->GetLocalRow();
236               colMerge += RecPoint2->GetLocalCol();
237               RecPoint2->SetEnergy(-1);     // Mark merged cluster
238             }
239             rowMerge /= (iSave + 1);
240             colMerge /= (iSave + 1);
241           }
242
243           Float_t smear[3];
244
245           // The position smearing in row-direction (uniform over pad width)            
246           Int_t row = (Int_t) ((rowMerge - row0) / rowPadSize);
247           smear[0]  = (row + gRandom->Rndm()) * rowPadSize + row0;
248
249           // The position smearing in rphi-direction (Gaussian)
250           smear[1] = 0;
251           do
252             smear[1] = gRandom->Gaus(colMerge,fRphiSigma);
253           while ((smear[1] < col0                        ) ||
254                  (smear[1] > col0 + nColMax * colPadSize));
255
256           // Time direction stays unchanged
257           smear[2] = time1;
258          
259           // Transform into local coordinates
260           smear[0] = (Int_t) ((smear[0] -  row0) /  rowPadSize);
261           smear[1] = (Int_t) ((smear[1] -  col0) /  colPadSize);
262           smear[2] = (Int_t) ((smear[2] - time0) / timeBinSize);
263
264           // Add the smeared cluster to the output array 
265           Int_t detector  = RecPoint1->GetDetector();
266           Int_t digits[3] = {0};
267           TRD->AddRecPoint(smear,digits,detector,0.0);
268
269         }
270
271         // Clear the temporary cluster-array and delete the cluster
272         Chamber->Delete();
273
274       }
275     }
276   }
277
278   printf("AliTRDclusterizerV0::MakeCluster -- ");
279   printf("Found %d points.\n",TRD->RecPoints()->GetEntries());
280   printf("AliTRDclusterizerV0::MakeCluster -- ");
281   printf("Fill the cluster tree.\n");
282   ClusterTree->Fill();
283
284   return kTRUE;
285
286 }