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