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