Do not use Atan, removed from ROOT too
[u/mrichter/AliRoot.git] / TRD / AliTRDv0.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.12  1999/11/02 16:35:56  fca
19 New version of TRD introduced
20
21 Revision 1.11  1999/11/01 20:41:51  fca
22 Added protections against using the wrong version of FRAME
23
24 Revision 1.10  1999/09/29 09:24:35  fca
25 Introduction of the Copyright and cvs Log
26
27 */
28
29 ///////////////////////////////////////////////////////////////////////////////
30 //                                                                           //
31 //  Transition Radiation Detector version 0 -- fast simulator                //
32 //                                                                           //
33 //Begin_Html
34 /*
35 <img src="picts/AliTRDfullClass.gif">
36 */
37 //End_Html
38 //                                                                           //
39 //                                                                           //
40 ///////////////////////////////////////////////////////////////////////////////
41
42 #include <TMath.h>
43 #include <TRandom.h>
44 #include <TVector.h>
45
46 #include "AliTRDv0.h"
47 #include "AliRun.h"
48 #include "AliMC.h"
49 #include "AliConst.h"
50   
51 ClassImp(AliTRDv0)
52
53 //_____________________________________________________________________________
54 AliTRDv0::AliTRDv0(const char *name, const char *title) 
55          :AliTRD(name, title) 
56 {
57   //
58   // Standard constructor for Transition Radiation Detector version 0
59   //
60
61   fIdSens     = 0;
62   fHitsOn     = 0;
63
64   fIdChamber1 = 0;
65   fIdChamber2 = 0;
66   fIdChamber3 = 0;
67
68   fRphiSigma  = 0;
69   fRphiDist   = 0;
70
71 }
72  
73 //_____________________________________________________________________________
74 void AliTRDv0::CreateGeometry()
75 {
76   //
77   // Create the GEANT geometry for the Transition Radiation Detector - Version 0
78   // This version covers the full azimuth. 
79   //
80
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;
84
85   // Define the chambers
86   AliTRD::CreateGeometry();
87
88 }
89
90 //_____________________________________________________________________________
91 void AliTRDv0::CreateMaterials()
92 {
93   //
94   // Create materials for the Transition Radiation Detector
95   //
96
97   AliTRD::CreateMaterials();
98
99 }
100
101 //_____________________________________________________________________________
102 void 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
287 //_____________________________________________________________________________
288 void AliTRDv0::Init() 
289 {
290   //
291   // Initialise Transition Radiation Detector after geometry is built
292   //
293
294   AliTRD::Init();
295
296   // Identifier of the sensitive volume (amplification region)
297   fIdSens     = gMC->VolId("UL06");
298
299   // Identifier of the TRD-driftchambers
300   fIdChamber1 = gMC->VolId("UCIO");
301   fIdChamber2 = gMC->VolId("UCIM");
302   fIdChamber3 = gMC->VolId("UCII");
303
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   
315 }
316
317 //_____________________________________________________________________________
318 void AliTRDv0::StepManager()
319 {
320   //
321   // Procedure called at every step in the TRD
322   // Fast simulator. If switched on, a hit is produced when a track
323   // crosses the border between amplification region and pad plane.
324   //
325
326   Int_t   vol[3]; 
327   Int_t   iIdSens, icSens; 
328   Int_t   iIdChamber, icChamber;
329
330   Float_t hits[4];
331
332   TLorentzVector p;
333   TClonesArray  &lhits = *fHits;
334
335   // Writing out hits enabled?
336   if (!(fHitsOn)) return;
337
338   // Use only charged tracks and count them only once per volume
339   if (gMC->TrackCharge()    && 
340       gMC->IsTrackExiting()) {
341     
342     // Check on sensitive volume
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
351       // The sector number
352       Float_t phi = hits[1] != 0 ? kRaddeg*TMath::ATan2(hits[0],hits[1]) : (hits[0] > 0 ? 180. : 0.);
353       vol[0] = ((Int_t) (phi / 20)) + 1;
354
355       // The chamber number 
356       //   1: outer left
357       //   2: middle left
358       //   3: inner
359       //   4: middle right
360       //   5: outer right
361       iIdChamber = gMC->CurrentVolOffID(1,icChamber);
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)       
367         vol[1] = 3;
368
369       // The plane number
370       vol[2] = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6;
371
372       new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
373
374     }
375
376   }  
377
378 }