TRDrec: installing headers
[u/mrichter/AliRoot.git] / TRD / TRDsim / AliTRDv1.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 //  Transition Radiation Detector version 1 -- slow simulator             //
21 //                                                                        //
22 ////////////////////////////////////////////////////////////////////////////
23
24 #include <TLorentzVector.h>
25 #include <TMath.h>
26 #include <TRandom.h>
27 #include <TVirtualMC.h>
28 #include <TGeoManager.h>
29 #include <TGeoMatrix.h>
30 #include <TGeoPhysicalNode.h>
31
32 #include "AliTrackReference.h"
33 #include "AliMC.h"
34 #include "AliRun.h"
35 #include "AliGeomManager.h"
36
37 #include "AliTRDgeometry.h"
38 #include "AliTRDCommonParam.h"
39 #include "AliTRDsimTR.h"
40 #include "AliTRDv1.h"
41
42 ClassImp(AliTRDv1)
43  
44 //_____________________________________________________________________________
45 AliTRDv1::AliTRDv1()
46   :AliTRD()
47   ,fTRon(kTRUE)
48   ,fTR(NULL)
49   ,fStepSize(0)
50   ,fWion(0)
51 {
52   //
53   // Default constructor
54   //
55
56 }
57
58 //_____________________________________________________________________________
59 AliTRDv1::AliTRDv1(const char *name, const char *title) 
60   :AliTRD(name,title) 
61   ,fTRon(kTRUE)
62   ,fTR(NULL)
63   ,fStepSize(0.1)
64   ,fWion(0)
65 {
66   //
67   // Standard constructor for Transition Radiation Detector version 1
68   //
69
70   SetBufferSize(128000);
71
72   if      (AliTRDCommonParam::Instance()->IsXenon()) {
73     fWion = 23.53; // Ionization energy XeCO2 (85/15)
74   }
75   else if (AliTRDCommonParam::Instance()->IsArgon()) {
76     fWion = 27.21; // Ionization energy ArCO2 (82/18)
77   }
78   else {
79     AliFatal("Wrong gas mixture");
80     exit(1);
81   }
82
83 }
84
85 //_____________________________________________________________________________
86 AliTRDv1::~AliTRDv1()
87 {
88   //
89   // AliTRDv1 destructor
90   //
91
92   if (fTR) {
93     delete fTR;
94     fTR     = 0;
95   }
96
97 }
98  
99 //_____________________________________________________________________________
100 void AliTRDv1::AddAlignableVolumes() const
101 {
102   //
103   // Create entries for alignable volumes associating the symbolic volume
104   // name with the corresponding volume path. Needs to be syncronized with
105   // eventual changes in the geometry.
106   //
107
108   TString volPath;
109   TString symName;
110
111   TString vpStr   = "ALIC_1/B077_1/BSEGMO";
112   TString vpApp1  = "_1/BTRD";
113   TString vpApp2  = "_1";
114   TString vpApp3a = "/UTR1_1/UTS1_1/UTI1_1/UT";
115   TString vpApp3b = "/UTR2_1/UTS2_1/UTI2_1/UT";
116   TString vpApp3c = "/UTR3_1/UTS3_1/UTI3_1/UT";
117   TString vpApp3d = "/UTR4_1/UTS4_1/UTI4_1/UT";
118
119   TString snStr   = "TRD/sm";
120   TString snApp1  = "/st";
121   TString snApp2  = "/pl";
122
123   //
124   // The super modules
125   // The symbolic names are: TRD/sm00
126   //                           ...
127   //                         TRD/sm17
128   //
129   for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
130
131     volPath  = vpStr;
132     volPath += isector;
133     volPath += vpApp1;
134     volPath += isector;
135     volPath += vpApp2;
136
137     symName  = snStr;
138     symName += Form("%02d",isector);
139
140     gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
141
142   }
143
144   //
145   // The readout chambers
146   // The symbolic names are: TRD/sm00/st0/pl0
147   //                           ...
148   //                         TRD/sm17/st4/pl5
149   //
150   AliGeomManager::ELayerID idTRD1 = AliGeomManager::kTRD1;
151   Int_t layer, modUID;
152   
153   for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
154
155     if (fGeometry->GetSMstatus(isector) == 0) continue;
156
157     for (Int_t istack = 0; istack < AliTRDgeometry::Nstack(); istack++) {
158       for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) {
159
160         layer = idTRD1 + ilayer;
161         modUID = AliGeomManager::LayerToVolUIDSafe(layer,isector*5+istack);
162
163         Int_t idet = AliTRDgeometry::GetDetectorSec(ilayer,istack);
164
165         volPath  = vpStr;
166         volPath += isector;
167         volPath += vpApp1;
168         volPath += isector;
169         volPath += vpApp2;
170         switch (isector) {
171         case 17:
172           if ((istack == 4) && (ilayer == 4)) {
173             continue;
174           }
175           volPath += vpApp3d;
176           break;
177         case 13:
178         case 14:
179         case 15:
180           if (istack == 2) {
181             continue;
182           }
183           volPath += vpApp3c;
184           break;
185         case 11:
186         case 12:
187           volPath += vpApp3b;
188           break;
189         default:
190           volPath += vpApp3a;
191         };
192         volPath += Form("%02d",idet);
193         volPath += vpApp2;
194
195         symName  = snStr;
196         symName += Form("%02d",isector);
197         symName += snApp1;
198         symName += istack;
199         symName += snApp2;
200         symName += ilayer;
201
202         TGeoPNEntry *alignableEntry = 
203           gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data(),modUID);
204
205         // Add the tracking to local matrix following the TPC example
206         if (alignableEntry) {
207           TGeoHMatrix *globMatrix = alignableEntry->GetGlobalOrig();
208           Double_t sectorAngle = 20.0 * (isector % 18) + 10.0;
209           TGeoHMatrix *t2lMatrix  = new TGeoHMatrix();
210           t2lMatrix->RotateZ(sectorAngle);
211           t2lMatrix->MultiplyLeft(&(globMatrix->Inverse()));
212           alignableEntry->SetMatrix(t2lMatrix);
213         }
214         else {
215           AliError(Form("Alignable entry %s is not valid!",symName.Data()));
216         }
217
218       }
219     }
220   }
221
222 }
223
224 //_____________________________________________________________________________
225 void AliTRDv1::CreateGeometry()
226 {
227   //
228   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
229   // This version covers the full azimuth. 
230   //
231
232   // Check that FRAME is there otherwise we have no place where to put the TRD
233   AliModule* frame = gAlice->GetModule("FRAME");
234   if (!frame) {
235     AliError("TRD needs FRAME to be present\n");
236     return;
237   }
238
239   // Define the chambers
240   AliTRD::CreateGeometry();
241
242 }
243
244 //_____________________________________________________________________________
245 void AliTRDv1::CreateMaterials()
246 {
247   //
248   // Create materials for the Transition Radiation Detector version 1
249   //
250
251   AliTRD::CreateMaterials();
252
253 }
254
255 //_____________________________________________________________________________
256 void AliTRDv1::CreateTRhit(Int_t det)
257 {
258   //
259   // Creates an electron cluster from a TR photon.
260   // The photon is assumed to be created a the end of the radiator. The 
261   // distance after which it deposits its energy takes into account the 
262   // absorbtion of the entrance window and of the gas mixture in drift
263   // volume.
264   //
265
266   // Maximum number of TR photons per track
267   const Int_t   kNTR         = 50;
268
269   TLorentzVector mom;
270   TLorentzVector pos;
271
272   Float_t eTR[kNTR];
273   Int_t   nTR;
274
275   // Create TR photons
276   TVirtualMC::GetMC()->TrackMomentum(mom);
277   Float_t pTot = mom.Rho();
278   fTR->CreatePhotons(11,pTot,nTR,eTR);
279   if (nTR > kNTR) {
280     AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
281   }
282
283   // Loop through the TR photons
284   for (Int_t iTR = 0; iTR < nTR; iTR++) {
285
286     Float_t energyMeV = eTR[iTR] * 0.001;
287     Float_t energyeV  = eTR[iTR] * 1000.0;
288     Float_t absLength = 0.0;
289     Float_t sigma     = 0.0;
290
291     // Take the absorbtion in the entrance window into account
292     Double_t muMy = fTR->GetMuMy(energyMeV);
293     sigma         = muMy * fFoilDensity;
294     if (sigma > 0.0) {
295       absLength = gRandom->Exp(1.0/sigma);
296       if (absLength < AliTRDgeometry::MyThick()) {
297         continue;
298       }
299     }
300     else {
301       continue;
302     }
303
304     // The absorbtion cross sections in the drift gas
305     // Gas-mixture (Xe/CO2)
306     Double_t muNo = 0.0;
307     if      (AliTRDCommonParam::Instance()->IsXenon()) {
308       muNo = fTR->GetMuXe(energyMeV);
309     }
310     else if (AliTRDCommonParam::Instance()->IsArgon()) {
311       muNo = fTR->GetMuAr(energyMeV);
312     }
313     Double_t muCO = fTR->GetMuCO(energyMeV);
314     sigma = (fGasNobleFraction * muNo + (1.0 - fGasNobleFraction) * muCO) 
315           * fGasDensity 
316           * fTR->GetTemp();
317
318     // The distance after which the energy of the TR photon
319     // is deposited.
320     if (sigma > 0.0) {
321       absLength = gRandom->Exp(1.0/sigma);
322       if (absLength > (AliTRDgeometry::DrThick()
323                      + AliTRDgeometry::AmThick())) {
324         continue;
325       }
326     }
327     else {
328       continue;
329     }
330
331     // The position of the absorbtion
332     Float_t posHit[3];
333     TVirtualMC::GetMC()->TrackPosition(pos);
334     posHit[0] = pos[0] + mom[0] / pTot * absLength;
335     posHit[1] = pos[1] + mom[1] / pTot * absLength;
336     posHit[2] = pos[2] + mom[2] / pTot * absLength;
337
338     // Create the charge 
339     Int_t q = ((Int_t) (energyeV / fWion));
340
341     // Add the hit to the array. TR photon hits are marked 
342     // by negative charge
343     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
344           ,det
345           ,posHit
346           ,-q
347           ,TVirtualMC::GetMC()->TrackTime()*1.0e06
348           ,kTRUE);
349
350   }
351
352 }
353
354 //_____________________________________________________________________________
355 void AliTRDv1::Init() 
356 {
357   //
358   // Initialise Transition Radiation Detector after geometry has been built.
359   //
360
361   AliTRD::Init();
362
363   AliDebug(1,"Slow simulator\n");
364
365   // Switch on TR simulation as default
366   if (!fTRon) {
367     AliInfo("TR simulation off");
368   }
369   else {
370     fTR = new AliTRDsimTR();
371   }
372
373   AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
374
375 }
376
377 //_____________________________________________________________________________
378 void AliTRDv1::StepManager()
379 {
380   //
381   // Slow simulator. Every charged track produces electron cluster as hits 
382   // along its path across the drift volume. The step size is fixed in
383   // this version of the step manager.
384   //
385   // Works for Xe/CO2 as well as Ar/CO2
386   //
387
388   // PDG code electron
389   const Int_t   kPdgElectron = 11;
390
391   Int_t    layer  = 0;
392   Int_t    stack  = 0;
393   Int_t    sector = 0;
394   Int_t    det    = 0;
395   Int_t    qTot;
396
397   Float_t  hits[3];
398   Double_t eDep;
399
400   Bool_t   drRegion = kFALSE;
401   Bool_t   amRegion = kFALSE;
402
403   TString  cIdPath;
404   Char_t   cIdSector[3];
405            cIdSector[2]  = 0;
406
407   TString  cIdCurrent;
408   TString  cIdSensDr = "J";
409   TString  cIdSensAm = "K";
410   Char_t   cIdChamber[3];
411            cIdChamber[2] = 0;
412
413   TLorentzVector pos;
414   TLorentzVector mom;
415
416   const Int_t    kNlayer      = AliTRDgeometry::Nlayer();
417   const Int_t    kNstack      = AliTRDgeometry::Nstack();
418   const Int_t    kNdetsec     = kNlayer * kNstack;
419
420   const Double_t kBig         = 1.0e+12;
421   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
422
423   // Set the maximum step size to a very large number for all 
424   // neutral particles and those outside the driftvolume
425   if (!fPrimaryIonisation) TVirtualMC::GetMC()->SetMaxStep(kBig); 
426
427   // If not charged track or already stopped or disappeared, just return.
428   if ((!TVirtualMC::GetMC()->TrackCharge()) || 
429         TVirtualMC::GetMC()->IsTrackDisappeared()) {
430     return;
431   }
432
433   // Inside a sensitive volume?
434   cIdCurrent = TVirtualMC::GetMC()->CurrentVolName();
435
436   if (cIdSensDr == cIdCurrent[1]) {
437     drRegion = kTRUE;
438   }
439   if (cIdSensAm == cIdCurrent[1]) {
440     amRegion = kTRUE;
441   }
442
443   if ((!drRegion) && 
444       (!amRegion)) {
445     return;
446   }
447
448   // The hit coordinates and charge
449   TVirtualMC::GetMC()->TrackPosition(pos);
450   hits[0] = pos[0];
451   hits[1] = pos[1];
452   hits[2] = pos[2];
453
454   // The sector number (0 - 17), according to standard coordinate system
455   cIdPath      = gGeoManager->GetPath();
456   cIdSector[0] = cIdPath[21];
457   cIdSector[1] = cIdPath[22];
458   sector = atoi(cIdSector);
459
460   // The plane and chamber number
461   cIdChamber[0]   = cIdCurrent[2];
462   cIdChamber[1]   = cIdCurrent[3];
463   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
464   stack = ((Int_t) idChamber / kNlayer);
465   layer = ((Int_t) idChamber % kNlayer);
466
467   // The detector number
468   det = fGeometry->GetDetector(layer,stack,sector);
469
470   // 0: InFlight 1:Entering 2:Exiting
471   Int_t trkStat = 0;
472
473   // Special hits only in the drift region
474   if      ((drRegion) &&
475            (TVirtualMC::GetMC()->IsTrackEntering())) {
476
477     // Create a track reference at the entrance of each
478     // chamber that contains the momentum components of the particle
479     TVirtualMC::GetMC()->TrackMomentum(mom);
480     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
481     trkStat = 1;
482
483     // Create the hits from TR photons if electron/positron is
484     // entering the drift volume
485     if ((fTR)   &&
486         (fTRon) &&
487         (TMath::Abs(TVirtualMC::GetMC()->TrackPid()) == kPdgElectron)) {
488       CreateTRhit(det);
489     }
490
491   }
492   else if ((amRegion) && 
493            (TVirtualMC::GetMC()->IsTrackExiting())) {
494
495     // Create a track reference at the exit of each
496     // chamber that contains the momentum components of the particle
497     TVirtualMC::GetMC()->TrackMomentum(mom);
498     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
499     trkStat = 2;
500
501   }
502   
503   // Calculate the charge according to GEANT Edep
504   // Create a new dEdx hit
505   eDep = TMath::Max(TVirtualMC::GetMC()->Edep(),0.0) * 1.0e+09;
506   qTot = (Int_t) (eDep / fWion);
507   if ((qTot) ||
508       (trkStat)) {
509     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
510           ,det
511           ,hits
512           ,qTot
513           ,TVirtualMC::GetMC()->TrackTime()*1.0e06
514           ,drRegion);
515   }
516
517   // Set Maximum Step Size
518   // Produce only one hit if Ekin is below cutoff
519   if ((TVirtualMC::GetMC()->Etot() - TVirtualMC::GetMC()->TrackMass()) < kEkinMinStep) {
520     return;
521   }
522   if (!fPrimaryIonisation) TVirtualMC::GetMC()->SetMaxStep(fStepSize);
523
524 }