MC-dependent part of AliRun extracted in AliMC (F.Carminati)
[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 //Begin_Html
23 /*
24 <img src="picts/AliTRDfullClass.gif">
25 */
26 //End_Html
27 //                                                                           //
28 //                                                                           //
29 ///////////////////////////////////////////////////////////////////////////////
30
31 #include <stdlib.h> 
32
33 #include <TF1.h>
34 #include <TLorentzVector.h>
35 #include <TMath.h>
36 #include <TRandom.h>
37 #include <TVector.h>
38 #include <TVirtualMC.h>
39
40 #include "AliConst.h"
41 #include "AliRun.h"
42 #include "AliTRDgeometry.h"
43 #include "AliTRDhit.h"
44 #include "AliTRDmatrix.h"
45 #include "AliTRDsim.h"
46 #include "AliTRDv1.h"
47 #include "AliMC.h"
48
49 ClassImp(AliTRDv1)
50  
51 //_____________________________________________________________________________
52 AliTRDv1::AliTRDv1():AliTRD()
53 {
54   //
55   // Default constructor
56   //
57
58   fSensSelect      =  0;
59   fSensPlane       = -1;
60   fSensChamber     = -1;
61   fSensSector      = -1;
62   fSensSectorRange =  0;
63
64   fDeltaE          = NULL;
65   fTR              = NULL;
66
67 }
68
69 //_____________________________________________________________________________
70 AliTRDv1::AliTRDv1(const char *name, const char *title) 
71          :AliTRD(name, title) 
72 {
73   //
74   // Standard constructor for Transition Radiation Detector version 1
75   //
76
77   fSensSelect      =  0;
78   fSensPlane       = -1;
79   fSensChamber     = -1;
80   fSensSector      = -1;
81   fSensSectorRange =  0;
82
83   fDeltaE          = NULL;
84   fTR              = NULL;
85
86   SetBufferSize(128000);
87
88 }
89
90 //_____________________________________________________________________________
91 AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
92 {
93   //
94   // Copy constructor
95   //
96
97   ((AliTRDv1 &) trd).Copy(*this);
98
99 }
100
101 //_____________________________________________________________________________
102 AliTRDv1::~AliTRDv1()
103 {
104   //
105   // AliTRDv1 destructor
106   //
107
108   if (fDeltaE) delete fDeltaE;
109   if (fTR)     delete fTR;
110
111 }
112  
113 //_____________________________________________________________________________
114 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
115 {
116   //
117   // Assignment operator
118   //
119
120   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
121   return *this;
122
123 }
124  
125 //_____________________________________________________________________________
126 void AliTRDv1::Copy(TObject &trd)
127 {
128   //
129   // Copy function
130   //
131
132   ((AliTRDv1 &) trd).fSensSelect      = fSensSelect;
133   ((AliTRDv1 &) trd).fSensPlane       = fSensPlane;
134   ((AliTRDv1 &) trd).fSensChamber     = fSensChamber;
135   ((AliTRDv1 &) trd).fSensSector      = fSensSector;
136   ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
137
138   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
139   fTR->Copy(*((AliTRDv1 &) trd).fTR);
140
141 }
142
143 //_____________________________________________________________________________
144 void AliTRDv1::CreateGeometry()
145 {
146   //
147   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
148   // This version covers the full azimuth. 
149   //
150
151   // Check that FRAME is there otherwise we have no place where to put the TRD
152   AliModule* frame = gAlice->GetModule("FRAME");
153   if (!frame) return;
154
155   // Define the chambers
156   AliTRD::CreateGeometry();
157
158 }
159
160 //_____________________________________________________________________________
161 void AliTRDv1::CreateMaterials()
162 {
163   //
164   // Create materials for the Transition Radiation Detector version 1
165   //
166
167   AliTRD::CreateMaterials();
168
169 }
170
171 //_____________________________________________________________________________
172 void AliTRDv1::CreateTRhit(Int_t det)
173 {
174   //
175   // Creates an electron cluster from a TR photon.
176   // The photon is assumed to be created a the end of the radiator. The 
177   // distance after which it deposits its energy takes into account the 
178   // absorbtion of the entrance window and of the gas mixture in drift
179   // volume.
180   //
181
182   // PDG code electron
183   const Int_t   kPdgElectron = 11;
184
185   // Ionization energy
186   const Float_t kWion        = 22.04;
187
188   // Maximum number of TR photons per track
189   const Int_t   kNTR         = 50;
190
191   TLorentzVector mom, pos;
192
193   // Create TR at the entrance of the chamber
194   if (gMC->IsTrackEntering()) {
195
196     // Create TR only for electrons 
197     Int_t iPdg = gMC->TrackPid();
198     if (TMath::Abs(iPdg) != kPdgElectron) return;
199
200     Float_t eTR[kNTR];
201     Int_t   nTR;
202
203     // Create TR photons
204     gMC->TrackMomentum(mom);
205     Float_t pTot = mom.Rho();
206     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
207     if (nTR > kNTR) {
208       printf("AliTRDv1::CreateTRhit -- ");
209       printf("Boundary error: nTR = %d, kNTR = %d\n",nTR,kNTR);
210       exit(1);
211     }
212
213     // Loop through the TR photons
214     for (Int_t iTR = 0; iTR < nTR; iTR++) {
215
216       Float_t energyMeV = eTR[iTR] * 0.001;
217       Float_t energyeV  = eTR[iTR] * 1000.0;
218       Float_t absLength = 0;
219       Float_t sigma     = 0;
220
221       // Take the absorbtion in the entrance window into account
222       Double_t muMy = fTR->GetMuMy(energyMeV);
223       sigma = muMy * fFoilDensity;
224       if (sigma > 0.0) {
225         absLength = gRandom->Exp(1.0/sigma);
226         if (absLength < AliTRDgeometry::MyThick()) continue;
227       }
228       else {
229         continue;
230       }
231
232       // The absorbtion cross sections in the drift gas
233       if (fGasMix == 1) {
234         // Gas-mixture (Xe/CO2)
235         Double_t muXe = fTR->GetMuXe(energyMeV);
236         Double_t muCO = fTR->GetMuCO(energyMeV);
237         sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
238       }
239       else {
240         // Gas-mixture (Xe/Isobutane) 
241         Double_t muXe = fTR->GetMuXe(energyMeV);
242         Double_t muBu = fTR->GetMuBu(energyMeV);
243         sigma = (0.97 * muXe + 0.03 * muBu) * fGasDensity * fTR->GetTemp();
244       }
245
246       // The distance after which the energy of the TR photon
247       // is deposited.
248       if (sigma > 0.0) {
249         absLength = gRandom->Exp(1.0/sigma);
250         if (absLength > AliTRDgeometry::DrThick()) continue;
251       }
252       else {
253         continue;
254       }
255
256       // The position of the absorbtion
257       Float_t posHit[3];
258       gMC->TrackPosition(pos);
259       posHit[0] = pos[0] + mom[0] / pTot * absLength;
260       posHit[1] = pos[1] + mom[1] / pTot * absLength;
261       posHit[2] = pos[2] + mom[2] / pTot * absLength;      
262
263       // Create the charge 
264       Int_t q = ((Int_t) (energyeV / kWion));
265
266       // Add the hit to the array. TR photon hits are marked 
267       // by negative charge
268       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE); 
269
270     }
271
272   }
273
274 }
275
276 //_____________________________________________________________________________
277 void AliTRDv1::Init() 
278 {
279   //
280   // Initialise Transition Radiation Detector after geometry has been built.
281   //
282
283   AliTRD::Init();
284
285   if(fDebug) printf("%s: Slow simulator\n",ClassName());
286   if (fSensSelect) {
287     if (fSensPlane   >= 0)
288       printf("          Only plane %d is sensitive\n",fSensPlane);
289     if (fSensChamber >= 0)   
290       printf("          Only chamber %d is sensitive\n",fSensChamber);
291     if (fSensSector  >= 0) {
292       Int_t sens1  = fSensSector;
293       Int_t sens2  = fSensSector + fSensSectorRange;
294             sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
295                    * AliTRDgeometry::Nsect();
296       printf("          Only sectors %d - %d are sensitive\n",sens1,sens2-1);
297     }
298   }
299   if (fTR) 
300     printf("%s: TR simulation on\n",ClassName());
301   else
302     printf("%s: TR simulation off\n",ClassName());
303   printf("\n");
304
305   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
306   const Float_t kPoti = 12.1;
307   // Maximum energy (50 keV);
308   const Float_t kEend = 50000.0;
309   // Ermilova distribution for the delta-ray spectrum
310   Float_t poti = TMath::Log(kPoti);
311   Float_t eEnd = TMath::Log(kEend);
312   fDeltaE = new TF1("deltae",Ermilova,poti,eEnd,0);
313
314   if(fDebug) {
315     printf("%s: ",ClassName());
316     for (Int_t i = 0; i < 80; i++) printf("*");
317     printf("\n");
318   }
319
320 }
321
322 //_____________________________________________________________________________
323 AliTRDsim *AliTRDv1::CreateTR()
324 {
325   //
326   // Enables the simulation of TR
327   //
328
329   fTR = new AliTRDsim();
330   return fTR;
331
332 }
333
334 //_____________________________________________________________________________
335 void AliTRDv1::SetSensPlane(Int_t iplane)
336 {
337   //
338   // Defines the hit-sensitive plane (0-5)
339   //
340
341   if ((iplane < 0) || (iplane > 5)) {
342     printf("Wrong input value: %d\n",iplane);
343     printf("Use standard setting\n");
344     fSensPlane  = -1;
345     fSensSelect =  0;
346     return;
347   }
348
349   fSensSelect = 1;
350   fSensPlane  = iplane;
351
352 }
353
354 //_____________________________________________________________________________
355 void AliTRDv1::SetSensChamber(Int_t ichamber)
356 {
357   //
358   // Defines the hit-sensitive chamber (0-4)
359   //
360
361   if ((ichamber < 0) || (ichamber > 4)) {
362     printf("Wrong input value: %d\n",ichamber);
363     printf("Use standard setting\n");
364     fSensChamber = -1;
365     fSensSelect  =  0;
366     return;
367   }
368
369   fSensSelect  = 1;
370   fSensChamber = ichamber;
371
372 }
373
374 //_____________________________________________________________________________
375 void AliTRDv1::SetSensSector(Int_t isector)
376 {
377   //
378   // Defines the hit-sensitive sector (0-17)
379   //
380
381   SetSensSector(isector,1);
382
383 }
384
385 //_____________________________________________________________________________
386 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
387 {
388   //
389   // Defines a range of hit-sensitive sectors. The range is defined by
390   // <isector> (0-17) as the starting point and <nsector> as the number 
391   // of sectors to be included.
392   //
393
394   if ((isector < 0) || (isector > 17)) {
395     printf("Wrong input value <isector>: %d\n",isector);
396     printf("Use standard setting\n");
397     fSensSector      = -1;
398     fSensSectorRange =  0;
399     fSensSelect      =  0;
400     return;
401   }
402
403   if ((nsector < 1) || (nsector > 18)) {
404     printf("Wrong input value <nsector>: %d\n",nsector);
405     printf("Use standard setting\n");
406     fSensSector      = -1;
407     fSensSectorRange =  0;
408     fSensSelect      =  0;
409     return;
410   }
411
412   fSensSelect      = 1;
413   fSensSector      = isector;
414   fSensSectorRange = nsector;
415
416 }
417
418 //_____________________________________________________________________________
419 void AliTRDv1::StepManager()
420 {
421   //
422   // Slow simulator. Every charged track produces electron cluster as hits 
423   // along its path across the drift volume. The step size is set acording
424   // to Bethe-Bloch. The energy distribution of the delta electrons follows
425   // a spectrum taken from Ermilova et al.
426   //
427
428   Int_t    pla = 0;
429   Int_t    cha = 0;
430   Int_t    sec = 0;
431   Int_t    det = 0;
432   Int_t    iPdg;
433   Int_t    qTot;
434
435   Float_t  hits[3];
436   Double_t  random[1];
437   Float_t  charge;
438   Float_t  aMass;
439
440   Double_t pTot = 0;
441   Double_t eDelta;
442   Double_t betaGamma, pp;
443   Double_t stepSize;
444
445   Bool_t   drRegion = kFALSE;
446   Bool_t   amRegion = kFALSE;
447
448   TString  cIdCurrent;
449   TString  cIdSensDr = "J";
450   TString  cIdSensAm = "K";
451   Char_t   cIdChamber[3];
452            cIdChamber[2] = 0;
453
454   TLorentzVector pos, mom;
455
456   const Int_t    kNplan       = AliTRDgeometry::Nplan();
457   const Int_t    kNcham       = AliTRDgeometry::Ncham();
458   const Int_t    kNdetsec     = kNplan * kNcham;
459
460   const Double_t kBig         = 1.0E+12;
461
462   // Ionization energy
463   const Float_t  kWion        = 22.04;
464   // Maximum momentum for e+ e- g 
465   const Float_t  kPTotMaxEl   = 0.002;
466   // Minimum energy for the step size adjustment
467   const Float_t  kEkinMinStep = 1.0e-5;
468   // Plateau value of the energy-loss for electron in xenon
469   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
470   //const Double_t kPlateau = 1.70;
471   // the averaged value (26/3/99)
472   const Float_t  kPlateau     = 1.55;
473   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
474   const Float_t  kPrim        = 48.0;
475   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
476   const Float_t  kPoti        = 12.1;
477
478   // PDG code electron
479   const Int_t    kPdgElectron = 11;
480
481   // Set the maximum step size to a very large number for all 
482   // neutral particles and those outside the driftvolume
483   gMC->SetMaxStep(kBig); 
484
485   // Use only charged tracks 
486   if (( gMC->TrackCharge()       ) &&
487       (!gMC->IsTrackStop()       ) && 
488       (!gMC->IsTrackDisappeared())) {
489
490     // Inside a sensitive volume?
491     drRegion = kFALSE;
492     amRegion = kFALSE;
493     cIdCurrent = gMC->CurrentVolName();
494     if (cIdSensDr == cIdCurrent[1]) {
495       drRegion = kTRUE;
496     }
497     if (cIdSensAm == cIdCurrent[1]) {
498       amRegion = kTRUE;
499     }
500     if (drRegion || amRegion) {
501
502       // The hit coordinates and charge
503       gMC->TrackPosition(pos);
504       hits[0] = pos[0];
505       hits[1] = pos[1];
506       hits[2] = pos[2];
507
508       // The sector number (0 - 17)
509       // The numbering goes clockwise and starts at y = 0
510       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
511       if (phi < 90.) 
512         phi = phi + 270.;
513       else
514         phi = phi -  90.;
515       sec = ((Int_t) (phi / 20));
516
517       // The plane and chamber number
518       cIdChamber[0] = cIdCurrent[2];
519       cIdChamber[1] = cIdCurrent[3];
520       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
521       cha = ((Int_t) idChamber / kNplan);
522       pla = ((Int_t) idChamber % kNplan);
523
524       // Check on selected volumes
525       Int_t addthishit = 1;
526       if (fSensSelect) {
527         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
528         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
529         if (fSensSector  >= 0) {
530           Int_t sens1  = fSensSector;
531           Int_t sens2  = fSensSector + fSensSectorRange;
532                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
533                        * AliTRDgeometry::Nsect();
534           if (sens1 < sens2) {
535             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
536           }
537           else {
538             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
539           }
540         }
541       }
542
543       // Add this hit
544       if (addthishit) {
545
546         // The detector number
547         det = fGeometry->GetDetector(pla,cha,sec);
548
549         // Special hits and TR photons only in the drift region
550         if (drRegion) {
551
552           // Create a track reference at the entrance and
553           // exit of each chamber that contain the 
554           // momentum components of the particle
555           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
556             gMC->TrackMomentum(mom);
557             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
558           }
559
560           // Create the hits from TR photons
561           if (fTR) CreateTRhit(det);
562
563         }
564
565         // Calculate the energy of the delta-electrons
566         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
567         eDelta = TMath::Max(eDelta,0.0);
568
569         // The number of secondary electrons created
570         qTot = ((Int_t) (eDelta / kWion) + 1);
571
572         // Create a new dEdx hit
573         if (drRegion) {
574           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,kTRUE);       
575         }
576         else {
577           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,kFALSE);      
578         }
579
580         // Calculate the maximum step size for the next tracking step
581         // Produce only one hit if Ekin is below cutoff 
582         aMass = gMC->TrackMass();
583         if ((gMC->Etot() - aMass) > kEkinMinStep) {
584
585           // The energy loss according to Bethe Bloch
586           iPdg  = TMath::Abs(gMC->TrackPid());
587           if ( (iPdg != kPdgElectron) ||
588               ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
589             gMC->TrackMomentum(mom);
590             pTot      = mom.Rho();
591             betaGamma = pTot / aMass;
592             pp        = kPrim * BetheBloch(betaGamma);
593             // Take charge > 1 into account
594             charge = gMC->TrackCharge();
595             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
596           }
597           // Electrons above 20 Mev/c are at the plateau
598           else {
599             pp = kPrim * kPlateau;
600           }
601       
602           if (pp > 0) {
603             do 
604             gMC->GetRandom()->RndmArray(1, random);
605             while ((random[0] == 1.) || (random[0] == 0.));
606             stepSize = - TMath::Log(random[0]) / pp; 
607             gMC->SetMaxStep(stepSize);
608           }
609
610         }
611
612       }
613
614     }
615
616   }
617
618 }
619
620 //_____________________________________________________________________________
621 Double_t AliTRDv1::BetheBloch(Double_t bg) 
622 {
623   //
624   // Parametrization of the Bethe-Bloch-curve
625   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
626   //
627
628   // This parameters have been adjusted to averaged values from GEANT
629   const Double_t kP1 = 7.17960e-02;
630   const Double_t kP2 = 8.54196;
631   const Double_t kP3 = 1.38065e-06;
632   const Double_t kP4 = 5.30972;
633   const Double_t kP5 = 2.83798;
634
635   // This parameters have been adjusted to Xe-data found in:
636   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
637   //const Double_t kP1 = 0.76176E-1;
638   //const Double_t kP2 = 10.632;
639   //const Double_t kP3 = 3.17983E-6;
640   //const Double_t kP4 = 1.8631;
641   //const Double_t kP5 = 1.9479;
642
643   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
644   const Double_t kBgMin = 0.8;
645   const Double_t kBBMax = 6.83298;
646   //const Double_t kBgMin = 0.6;
647   //const Double_t kBBMax = 17.2809;
648   //const Double_t kBgMin = 0.4;
649   //const Double_t kBBMax = 82.0;
650
651   if (bg > kBgMin) {
652     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
653     Double_t aa = TMath::Power(yy,kP4);
654     Double_t bb = TMath::Power((1./bg),kP5);
655              bb = TMath::Log(kP3 + bb);
656     return ((kP2 - aa - bb)*kP1 / aa);
657   }
658   else {
659     return kBBMax;
660   }
661
662 }
663
664 //_____________________________________________________________________________
665 Double_t Ermilova(Double_t *x, Double_t *)
666 {
667   //
668   // Calculates the delta-ray energy distribution according to Ermilova.
669   // Logarithmic scale !
670   //
671
672   Double_t energy;
673   Double_t dpos;
674   Double_t dnde;
675
676   Int_t    pos1, pos2;
677
678   const Int_t kNv = 31;
679
680   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120  
681                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
682                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
683                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
684                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
685                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
686                      ,11.5129 };
687
688   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
689                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
690                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
691                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
692                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
693                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
694                      ,  0.004 };
695
696   energy = x[0];
697
698   // Find the position 
699   pos1 = pos2 = 0;
700   dpos = 0;
701   do {
702     dpos = energy - vxe[pos2++];
703   } 
704   while (dpos > 0);
705   pos2--; 
706   if (pos2 > kNv) pos2 = kNv - 1;
707   pos1 = pos2 - 1;
708
709   // Differentiate between the sampling points
710   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
711
712   return dnde;
713
714 }