Do not use Atan, removed from ROOT too
[u/mrichter/AliRoot.git] / TRD / AliTRDv0.cxx
CommitLineData
4c039060 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$
90f8d287 18Revision 1.12 1999/11/02 16:35:56 fca
19New version of TRD introduced
20
5c7f4665 21Revision 1.11 1999/11/01 20:41:51 fca
22Added protections against using the wrong version of FRAME
23
ab76897d 24Revision 1.10 1999/09/29 09:24:35 fca
25Introduction of the Copyright and cvs Log
26
4c039060 27*/
28
fe4da5cc 29///////////////////////////////////////////////////////////////////////////////
30// //
5c7f4665 31// Transition Radiation Detector version 0 -- fast simulator //
fe4da5cc 32// //
33//Begin_Html
34/*
5c7f4665 35<img src="picts/AliTRDfullClass.gif">
fe4da5cc 36*/
37//End_Html
38// //
39// //
40///////////////////////////////////////////////////////////////////////////////
41
42#include <TMath.h>
43#include <TRandom.h>
44#include <TVector.h>
fe4da5cc 45
fe4da5cc 46#include "AliTRDv0.h"
47#include "AliRun.h"
48#include "AliMC.h"
49#include "AliConst.h"
50
51ClassImp(AliTRDv0)
52
53//_____________________________________________________________________________
54AliTRDv0::AliTRDv0(const char *name, const char *title)
55 :AliTRD(name, title)
56{
57 //
58 // Standard constructor for Transition Radiation Detector version 0
59 //
82bbf98a 60
61 fIdSens = 0;
62 fHitsOn = 0;
63
82bbf98a 64 fIdChamber1 = 0;
65 fIdChamber2 = 0;
66 fIdChamber3 = 0;
67
5c7f4665 68 fRphiSigma = 0;
69 fRphiDist = 0;
70
fe4da5cc 71}
72
73//_____________________________________________________________________________
74void AliTRDv0::CreateGeometry()
75{
76 //
82bbf98a 77 // Create the GEANT geometry for the Transition Radiation Detector - Version 0
78 // This version covers the full azimuth.
d3f347ff 79 //
d3f347ff 80
82bbf98a 81 // Check that FRAME is there otherwise we have no place where to put the TRD
82 AliModule* FRAME = gAlice->GetModule("FRAME");
83 if (!FRAME) return;
fe4da5cc 84
82bbf98a 85 // Define the chambers
86 AliTRD::CreateGeometry();
fe4da5cc 87
fe4da5cc 88}
89
90//_____________________________________________________________________________
91void AliTRDv0::CreateMaterials()
92{
93 //
94 // Create materials for the Transition Radiation Detector
95 //
82bbf98a 96
fe4da5cc 97 AliTRD::CreateMaterials();
82bbf98a 98
fe4da5cc 99}
100
5c7f4665 101//_____________________________________________________________________________
102void AliTRDv0::Hits2Clusters()
103{
104 // A simple cluster generator. It takes the hits from the
105 // fast simulator (one hit per plane) and transforms them
106 // into cluster, by applying position smearing and merging
107 // of nearby cluster. The searing is done uniformly in z-direction
108 // over the length of a readout pad. In rphi-direction a Gaussian
109 // smearing is applied with a sigma given by fRphiSigma.
110 // Clusters are considered as overlapping when they are closer in
111 // rphi-direction than the value defined in fRphiDist.
112 // Use the macro fastClusterCreate.C to create the cluster.
113
114 printf("AliTRDv0::Hits2Clusters -- Start creating cluster\n");
115
116 Int_t nBytes = 0;
117
118 AliTRDhit *TRDhit;
119
120 // Get the pointer to the hit tree
121 TTree *HitTree = gAlice->TreeH();
122 // Get the pointer to the reconstruction tree
123 TTree *ClusterTree = gAlice->TreeD();
124
125 TObjArray *Chamber = new TObjArray();
126
127 // Get the number of entries in the hit tree
128 // (Number of primary particles creating a hit somewhere)
129 Int_t nTrack = (Int_t) HitTree->GetEntries();
130
131 // Loop through all the chambers
132 for (Int_t icham = 0; icham < kNcham; icham++) {
133 for (Int_t iplan = 0; iplan < kNplan; iplan++) {
134 for (Int_t isect = 0; isect < kNsect; isect++) {
135
136 // Loop through all entries in the tree
137 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
138
139 gAlice->ResetHits();
140 nBytes += HitTree->GetEvent(iTrack);
141
142 // Get the number of hits in the TRD created by this particle
143 Int_t nHit = fHits->GetEntriesFast();
144
145 // Loop through the TRD hits
146 for (Int_t iHit = 0; iHit < nHit; iHit++) {
147
148 if (!(TRDhit = (AliTRDhit *) fHits->UncheckedAt(iHit)))
149 continue;
150
151 Float_t x = TRDhit->fX;
152 Float_t y = TRDhit->fY;
153 Float_t z = TRDhit->fZ;
154 Int_t track = TRDhit->fTrack;
155 Int_t plane = TRDhit->fPlane;
156 Int_t sector = TRDhit->fSector;
157 Int_t chamber = TRDhit->fChamber;
158
159 if ((sector != isect+1) ||
160 (plane != iplan+1) ||
161 (chamber != icham+1))
162 continue;
163
164 // Rotate the sectors on top of each other
165 Float_t phi = 2.0 * kPI / (Float_t) kNsect
166 * ((Float_t) sector - 0.5);
167 Float_t xRot = -x * TMath::Cos(phi) + y * TMath::Sin(phi);
168 Float_t yRot = x * TMath::Sin(phi) + y * TMath::Cos(phi);
169 Float_t zRot = z;
170
171 // Add this cluster to the temporary cluster-array for this chamber
172 Int_t tracks[3];
173 tracks[0] = track;
174 Int_t clusters[5];
175 clusters[0] = sector;
176 clusters[1] = chamber;
177 clusters[2] = plane;
178 clusters[3] = 0;
179 clusters[4] = 0;
180 Float_t position[3];
181 position[0] = zRot;
182 position[1] = yRot;
183 position[2] = xRot;
184 AliTRDcluster *Cluster = new AliTRDcluster(tracks,clusters,position);
185 Chamber->Add(Cluster);
186
187 }
188
189 }
190
191 // Loop through the temporary cluster-array
192 for (Int_t iClus1 = 0; iClus1 < Chamber->GetEntries(); iClus1++) {
193
194 AliTRDcluster *Cluster1 = (AliTRDcluster *) Chamber->UncheckedAt(iClus1);
195 Float_t x1 = Cluster1->fX;
196 Float_t y1 = Cluster1->fY;
197 Float_t z1 = Cluster1->fZ;
198
199 if (!(z1)) continue; // Skip marked cluster
200
201 const Int_t nSave = 2;
202 Int_t idxSave[nSave];
203 Int_t iSave = 0;
204
205 Int_t tracks[3];
206 tracks[0] = Cluster1->fTracks[0];
207
208 // Check the other cluster to see, whether there are close ones
209 for (Int_t iClus2 = iClus1 + 1; iClus2 < Chamber->GetEntries(); iClus2++) {
210 AliTRDcluster *Cluster2 = (AliTRDcluster *) Chamber->UncheckedAt(iClus2);
211 Float_t x2 = Cluster2->fX;
212 Float_t y2 = Cluster2->fY;
213 if ((TMath::Abs(x1 - x2) < fRowPadSize) ||
214 (TMath::Abs(y1 - y2) < fRphiDist)) {
215 if (iSave == nSave) {
216 printf("AliTRDv0::Hits2Clusters -- Boundary error: iSave = %d, nSave = %d\n"
217 ,iSave,nSave);
218 }
219 else {
220 idxSave[iSave] = iClus2;
221 tracks[iSave+1] = Cluster2->fTracks[0];
222 }
223 iSave++;
224 }
225 }
226
227 // Merge close cluster
228 Float_t yMerge = y1;
229 Float_t xMerge = x1;
230 if (iSave) {
231 for (Int_t iMerge = 0; iMerge < iSave; iMerge++) {
232 AliTRDcluster *Cluster2 = (AliTRDcluster *) Chamber->UncheckedAt(idxSave[iMerge]);
233 xMerge += Cluster2->fX;
234 yMerge += Cluster2->fY;
235 Cluster2->fZ = 0; // Mark merged cluster
236 }
237 xMerge /= (iSave + 1);
238 yMerge /= (iSave + 1);
239 }
240
241 // The position smearing in z-direction (uniform over pad width)
242 Int_t row = (Int_t) ((xMerge - fRow0[iplan][icham][isect]) / fRowPadSize);
243 Float_t xSmear = (row + gRandom->Rndm()) * fRowPadSize
244 + fRow0[iplan][icham][isect];
245
246 // The position smearing in rphi-direction (Gaussian)
247 Float_t ySmear = 0;
248 do
249 ySmear = gRandom->Gaus(yMerge,fRphiSigma);
250 while ((ySmear < fCol0[iplan]) ||
251 (ySmear > fCol0[iplan] + fColMax[iplan] * fColPadSize));
252
253 // Time direction stays unchanged
254 Float_t zSmear = z1;
255
256 Int_t clusters[5];
257 clusters[0] = Cluster1->fSector;
258 clusters[1] = Cluster1->fChamber;
259 clusters[2] = Cluster1->fPlane;
260 clusters[3] = 0;
261 clusters[4] = 0;
262 Float_t position[3];
263 // Rotate the sectors back into their real position
264 Float_t phi = 2*kPI / kNsect * ((Float_t) Cluster1->fSector - 0.5);
265 position[0] = -zSmear * TMath::Cos(phi) + ySmear * TMath::Sin(phi);
266 position[1] = zSmear * TMath::Sin(phi) + ySmear * TMath::Cos(phi);
267 position[2] = xSmear;
268
269 // Add the smeared cluster to the output array
270 AddCluster(tracks,clusters,position);
271
272 }
273
274 // Clear the temporary cluster-array and delete the cluster
275 Chamber->Delete();
276
277 }
278 }
279 }
280
281 printf("AliTRDv0::Hits2Clusters -- Found %d cluster\n",fClusters->GetEntries());
282 printf("AliTRDv0::Hits2Clusters -- Fill the cluster tree\n");
283 ClusterTree->Fill();
284
285}
286
fe4da5cc 287//_____________________________________________________________________________
288void AliTRDv0::Init()
289{
290 //
291 // Initialise Transition Radiation Detector after geometry is built
292 //
82bbf98a 293
fe4da5cc 294 AliTRD::Init();
82bbf98a 295
82bbf98a 296 // Identifier of the sensitive volume (amplification region)
297 fIdSens = gMC->VolId("UL06");
298
82bbf98a 299 // Identifier of the TRD-driftchambers
300 fIdChamber1 = gMC->VolId("UCIO");
301 fIdChamber2 = gMC->VolId("UCIM");
302 fIdChamber3 = gMC->VolId("UCII");
303
5c7f4665 304 // Parameter for Hits2Cluster
305
306 // Position resolution in rphi-direction
307 fRphiSigma = 0.02;
308 // Minimum distance of non-overlapping cluster
309 fRphiDist = 1.0;
310
311 printf(" Fast simulator\n");
312 for (Int_t i = 0; i < 80; i++) printf("*");
313 printf("\n");
314
fe4da5cc 315}
316
317//_____________________________________________________________________________
318void AliTRDv0::StepManager()
319{
320 //
321 // Procedure called at every step in the TRD
82bbf98a 322 // Fast simulator. If switched on, a hit is produced when a track
323 // crosses the border between amplification region and pad plane.
fe4da5cc 324 //
325
82bbf98a 326 Int_t vol[3];
327 Int_t iIdSens, icSens;
82bbf98a 328 Int_t iIdChamber, icChamber;
329
82bbf98a 330 Float_t hits[4];
fe4da5cc 331
0a6d8768 332 TLorentzVector p;
82bbf98a 333 TClonesArray &lhits = *fHits;
fe4da5cc 334
82bbf98a 335 // Writing out hits enabled?
336 if (!(fHitsOn)) return;
fe4da5cc 337
fe4da5cc 338 // Use only charged tracks and count them only once per volume
82bbf98a 339 if (gMC->TrackCharge() &&
340 gMC->IsTrackExiting()) {
fe4da5cc 341
342 // Check on sensitive volume
82bbf98a 343 iIdSens = gMC->CurrentVolID(icSens);
344 if (iIdSens == fIdSens) {
345
346 gMC->TrackPosition(p);
347 for (Int_t i = 0; i < 3; i++) hits[i] = p[i];
348 // No charge created
349 hits[3] = 0;
350
fe4da5cc 351 // The sector number
90f8d287 352 Float_t phi = hits[1] != 0 ? kRaddeg*TMath::ATan2(hits[0],hits[1]) : (hits[0] > 0 ? 180. : 0.);
5c7f4665 353 vol[0] = ((Int_t) (phi / 20)) + 1;
82bbf98a 354
d3f347ff 355 // The chamber number
356 // 1: outer left
82bbf98a 357 // 2: middle left
d3f347ff 358 // 3: inner
82bbf98a 359 // 4: middle right
d3f347ff 360 // 5: outer right
5c7f4665 361 iIdChamber = gMC->CurrentVolOffID(1,icChamber);
82bbf98a 362 if (iIdChamber == fIdChamber1)
363 vol[1] = (hits[2] < 0 ? 1 : 5);
364 else if (iIdChamber == fIdChamber2)
365 vol[1] = (hits[2] < 0 ? 2 : 4);
366 else if (iIdChamber == fIdChamber3)
d3f347ff 367 vol[1] = 3;
82bbf98a 368
fe4da5cc 369 // The plane number
82bbf98a 370 vol[2] = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6;
371
372 new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
d3f347ff 373
fe4da5cc 374 }
d3f347ff 375
fe4da5cc 376 }
d3f347ff 377
fe4da5cc 378}