Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / CRT / AliGenCRT.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.1.2.1  2002/10/10 14:40:31  hristov
19 Updating VirtualMC to v3-09-02
20
21 Revision 1.1  2002/10/07 11:25:28  gamez
22 First version, generation of cosmic muons on the surface
23
24
25 */
26
27 /////////////////////////////////////////////////////////////////////////////
28 //
29 //  Contain parametrizations to generate atmospheric muons, and also
30 //  to generate single muons and muon bundles at surface level.
31 //
32 //Begin_Html
33 /*
34 <img src="picts/AliGenCRTClass.gif">
35 </pre>
36 <br clear=left>
37 <font size=+2 color=red>
38 <p>The responsible person for this module is
39 <a href="mailto:Enrique.Gamez.Flores@cern.ch">Enrique Gamez</a>.
40 </font>
41 <pre>
42 */
43 //End_Html
44 //
45 /////////////////////////////////////////////////////////////////////////////
46
47 #include <iostream.h>
48
49 #include "AliRun.h"
50 #include "AliConst.h"
51 #include "AliPDG.h"
52 #include "AliMCProcess.h"
53
54 #include "AliGenCRT.h"
55
56 ClassImp(AliGenCRT)
57
58 //_____________________________________________________________________________
59 AliGenCRT::AliGenCRT() : AliGenerator(-1)
60 {
61   //
62   // Default ctor.
63   //
64   fIpart = 0;
65
66   fXwidth=0;
67   fNx=1;
68   fZwidth=0;
69   fNz=1;
70   fMuonGrid = kFALSE;
71
72
73   // Set the origin above the vertex, on the surface.
74   fOrigin[0] = 0.;
75   fOrigin[1] = 0.;
76   fOrigin[2] = 0.;
77
78   fZenithMin = 0.; 
79   fZenithMax = 0.;
80
81   fAzimuthMin = 0.;
82   fAzimuthMax = 0.;
83
84   fPResolution = 1.;
85   fPRange = 0.;
86
87   fMomentumDist = 0;
88   fZenithDist = 0;
89
90   fPDist = 0;
91
92   fAp = 0;
93   fUnfoldedMomentumDist = 0;
94
95   fCRModeName = 0;
96 }
97
98 //_____________________________________________________________________________
99 AliGenCRT::AliGenCRT(Int_t npart) 
100   : AliGenerator(npart)
101 {
102   //
103   // Standard ctor.
104   //
105   fName = "CRT";
106   fTitle = "Cosmic Muons generator";
107
108   // Generate Muon- by default
109   fIpart = kMuonMinus;
110
111   fXwidth=0;
112   fNx=1;
113   fZwidth=0;
114   fNz=1;
115   fMuonGrid = kFALSE;
116
117   // Set the origin above the vertex, on the surface.
118   fOrigin[0] = 0.;
119   fOrigin[1] = AliCRTConstants::fgDepth; // At the surface by default.
120   fOrigin[2] = 0.;
121
122   fZenithMin = 0.; // Generate veritcals by default.
123   fZenithMax = 0.;
124
125   fAzimuthMin = 0.;
126   fAzimuthMax = 0.;
127
128   fPResolution = 1.; // 1 GeV by default.
129   fPRange = 0.;      // 0 GeV by default.
130
131   fMomentumDist = 0;
132   fZenithDist = 0;
133
134   fPDist = 0;
135
136   fAp = 0;
137   fUnfoldedMomentumDist = 0;
138
139   fCRModeName = 0;
140 }
141
142 //_____________________________________________________________________________
143 AliGenCRT::AliGenCRT(const AliGenCRT& gen) : AliGenerator()
144 {
145   //
146   // Copy ctor.
147   //
148   gen.Copy(*this);
149 }
150
151 //_____________________________________________________________________________
152 AliGenCRT& AliGenCRT::operator= (const AliGenCRT& gen)
153 {
154   //
155   // Asingment ctor.
156   //
157   gen.Copy(*this);
158   return *this;
159 }
160
161 //_____________________________________________________________________________
162 AliGenCRT::~AliGenCRT()
163 {
164   //
165   // Default dtor.
166   //
167   if ( fAp )           delete fAp;
168   if ( fMomentumDist ) delete fMomentumDist;
169   if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
170   if ( fZenithDist )   delete fZenithDist;
171   if ( fPDist )        delete fPDist;
172 }
173
174 //_____________________________________________________________________________
175 void AliGenCRT::Generate()
176 {
177   //
178   // Generate on one trigger
179   //
180
181   Float_t polar[3]= {0,0,0}; // Polarization parameters
182   //
183   Float_t origin[3];
184   Float_t p[3];
185   Int_t i, j, nt;
186   Float_t pmom, pt;
187   Float_t zenith, azimuth;
188   //
189   Float_t random[6];
190
191
192   // Check if you need to use a distribution function for the momentum
193   if ( fMomentumDist ) {
194     pmom = - this->GetMomentum(); // Always downwards.
195
196   } else {
197     pmom = -fPMin;
198
199   }
200   
201   zenith = this->GetZenithAngle(pmom);  // In degrees
202   
203   // Take the azimuth random.
204   Rndm(random, 2);
205   azimuth = (fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0]); // In degrees
206   
207   origin[0] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
208     TMath::Sin(azimuth*kDegrad);
209   origin[1] = AliCRTConstants::fgDepth;
210   origin[2] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
211     TMath::Cos(azimuth*kDegrad);
212   
213   if ( fVertexSmear == kPerEvent ) {
214     Rndm(random,6);
215     for (j=0;j<3;j++) {
216       if ( j == 1 ) {
217         // Don't smear the vertical position.
218         origin[j] = AliCRTConstants::fgDepth;
219       } else {
220         origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
221           TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
222       }
223     }
224   }
225
226
227   if ( fCRMode == kSingleMuons ) {
228
229     //
230     // Generate kinematics a la AliGenBox
231     //
232
233     for(i=0;i<fNpart;i++) {
234       Rndm(random,3);
235
236       if(TestBit(kMomentumRange)) {
237         pmom = -( fPMin + random[1]*(fPMax - fPMin) ); // always downwards.
238         pt=pmom*TMath::Sin(zenith*kDegrad);
239       } else {
240
241         pt= -(fPtMin+random[1]*(fPtMax-fPtMin)); // always downwards.
242         pmom=pt/TMath::Sin(zenith*kDegrad);
243       }
244       
245       p[0] = pt*TMath::Sin(azimuth*kDegrad);
246       p[1] = pmom*TMath::Cos(zenith*kDegrad);
247       p[2] = pt*TMath::Cos(azimuth*kDegrad);
248       
249       if(fVertexSmear==kPerTrack) {
250         Rndm(random,6);
251         for (j=0;j<3;j++) {
252           if ( j == 1 ) {
253             origin[j] = AliCRTConstants::fgDepth;
254           } else {
255             origin[j]=fOrigin[j]+fOsigma[j]*
256               TMath::Cos(2*random[2*j]*TMath::Pi())*
257               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
258           }
259         }
260       }
261
262     }
263     
264     // End of generation a la AliGenBox
265     
266   } else if ( fCRMode == kMuonBundle ) {
267
268     //
269     // Generate kinematics a la AliGenScan
270     //
271     Float_t dx,dz;
272     //
273     if (fNx > 0) {
274       dx=fXwidth/fNx;
275     } else {
276       dx=1e10;
277     }
278     
279     if (fNz > 0) {
280       dz=fZwidth/fNz;
281     } else {
282       dz=1e10;
283     }
284
285     for (Int_t ix=0; ix<fNx; ix++) {
286       for (Int_t iz=0; iz<fNz; iz++){
287         Rndm(random,6);
288         origin[0]+=ix*dx+2*(random[0]-0.5)*fOsigma[0];
289         origin[2]+=iz*dz+2*(random[1]-0.5)*fOsigma[2];       
290
291         pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downward.
292         p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
293         p[1] = TMath::Cos(zenith*kDegrad)*pmom;
294         p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
295         
296       }
297     }
298
299     // End of generation a la AliGenScan
300
301   } else if ( fCRMode == kMuonFlux ) {
302     
303     // Always generate the things downwards.
304     p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
305     p[1] = TMath::Cos(zenith*kDegrad)*pmom;
306     p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
307
308   } else {
309     Error("Generate", "Generation Mode unknown!\n");
310   }
311
312   // Put the track on the stack.
313   SetTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
314
315
316 }
317
318 //_____________________________________________________________________________
319 void AliGenCRT::Init()
320 {
321   //
322   // Initialize some internal methods.
323   //
324
325   // Determine some specific data members.
326   fPRange = TMath::Abs(fPMax-fPMin);
327
328   if ( fCRMode == kSingleMuons ) {
329     fCRModeName = new TString("Single Muons");
330     // Initialisation, check consistency of selected ranges
331     if(TestBit(kPtRange)&&TestBit(kMomentumRange)) 
332       Fatal("Init","You should not set the momentum range and the pt range!");
333     
334     if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange))) 
335       Fatal("Init","You should set either the momentum or the pt range!");
336     
337   } else if ( fCRMode == kMuonBundle ) {
338     fCRModeName = new TString("Muon Bundles");
339
340   } else if ( fCRMode == kMuonFlux ) {
341     fCRModeName = new TString("Muon Fluxes");
342     // Initialize the ditribution functions.
343     this->InitMomentumGeneration();
344     this->InitZenithalAngleGeneration();
345     
346   } else {
347     Fatal("Generate", "Generation Mode unknown!\n");
348
349   }
350
351 }
352
353 //____________________________________________________________________________
354 void AliGenCRT::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
355 {
356   //
357   // Define the grid
358   // This data shuold be used for Muon bundles generation.
359   //
360   fXwidth=xwidth;
361   fNx=nx;
362   fZwidth=zwidth;
363   fNz=nz;
364
365   // Print a message  about the use, if the Mode has not been set, or
366   // it has to a different Mode.
367   if ( fCRMode != kMuonBundle ) {
368     Warning("SetRange","You have been specified a grid to generate muon bundles, but seems that you haven't choose this generation mode, or you have already select a different one");
369     fMuonGrid = kTRUE;
370   }
371 }
372
373 //____________________________________________________________________________
374 void AliGenCRT::InitApWeightFactors()
375 {
376   //
377   // This factors will be  to correct the zenithal angle distribution
378   // acording the momentum
379
380   //
381   // Fill the array for the flux zenith angle dependence.
382   // at the index 0 of fAp[] will be the "factor" if we have a muon
383   // of 0 GeV.
384   Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
385   Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};
386
387   // Get the information from the momentum
388   Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
389   // Initialize the Array of floats to hold the a(p) factors.
390   fAp = new TArrayF(pEnd);
391   
392   Int_t index = 0;
393
394   for (Int_t i = 0; i < pEnd; i++ ) {
395     Float_t currentP = ((Float_t)i)*fPResolution;
396     if ( currentP < p[1] )                          index = 0;
397     else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
398     else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
399     else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
400     else if ( currentP >= p[4] )                    index = 4;
401
402     Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
403                  (p[index+1] - p[index]) + a[index];
404     fAp->AddAt(ap, i);
405   }
406
407 }
408
409 //___________________________________________________________________________
410 void AliGenCRT::InitMomentumGeneration()
411 {
412   //
413   // Initialize a funtion to generate the momentum randomly
414   // acording this function.
415   //
416
417   // Check if we nned to initialize the function
418   if ( fPMin != fPMax ) {
419
420     // Check also if the function have been defined yet.
421     if ( !fMomentumDist ) {
422
423       // If not, use this function
424       const char* y      = "log10(x)";
425       
426       const char* h1Coef = "[0]*( %s*%s*%s/2 - (5*%s*%s/2) + 3*%s )";
427       const char* h2Coef = "[1]*( (-2*%s*%s*%s/3) + (3*%s*%s) - 10*%s/3 + 1 )";
428       const char* h3Coef = "[2]*( %s*%s*%s/6 - %s*%s/2 + %s/3 )";
429       const char* s2Coef = "[3]*( %s*%s*%s/3 - 2*%s*%s + 11*%s/3 - 2 )";
430       
431       const char* h = "%s + %s + %s + %s";
432       const char* flux = "pow(10., %s)";
433       const char* normalizedFlux = "0.86*x*x*x*pow(10., %s)";
434       const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
435       
436       char buffer1[1024];
437       char buffer2[1024];
438       char buffer3[1024];
439       char buffer4[1024];
440       char buffer5[1024];
441       char buffer6[1024];
442       char buffer7[1024];
443
444       sprintf(buffer1, h1Coef, y, y, y, y, y, y);
445       sprintf(buffer2, h2Coef, y, y, y, y, y, y);
446       sprintf(buffer3, h3Coef, y, y, y, y, y, y);
447       sprintf(buffer4, s2Coef, y, y, y, y, y, y);
448       
449       sprintf(buffer5, h, buffer1, buffer2, buffer3, buffer4);
450       
451       sprintf(buffer6, flux, buffer5);
452       
453       fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
454       sprintf(buffer7, normalizedFlux, buffer5);
455       fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
456       for (Int_t i = 0; i < 4; i++ ) {
457         fMomentumDist->SetParName(i, paramNames[i]);
458         fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
459       }
460       
461       fMomentumDist->SetParameter(0, 0.133);
462       fMomentumDist->SetParameter(1, -2.521);
463       fMomentumDist->SetParameter(2, -5.78);
464       fMomentumDist->SetParameter(3, -2.11);
465
466       fUnfoldedMomentumDist->SetParameter(0, 0.133);
467       fUnfoldedMomentumDist->SetParameter(1, -2.521);
468       fUnfoldedMomentumDist->SetParameter(2, -5.78);
469       fUnfoldedMomentumDist->SetParameter(3, -2.11);
470       
471     }
472
473   }
474
475 }
476
477 //____________________________________________________________________________
478 void AliGenCRT::InitZenithalAngleGeneration()
479 {
480   //
481   // Initalize a distribution function for the zenith angle.
482   // This angle will be obtained randomly acording this function.
483   // The generated angles  will been in degrees.
484
485   // Check if we need to create the function.
486   if ( fZenithMin != fZenithMax ) {
487
488     // Check also if another function have been defined.
489     if ( !fZenithDist ) {
490       
491       // initialize the momentum dependent coefficients, a(p) 
492       InitApWeightFactors();
493       
494       // Define the standard function.
495       char* zenithalDisributionFunction = "1 + [0]*(1 - cos(x*3.14159265358979312/180))";
496
497       Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
498       fPDist = new TClonesArray("TF1", pEnd);
499       TClonesArray &angle = *fPDist;
500       for ( Int_t i = 0; i < pEnd; i++ ) {
501         // Fill the distribution
502         TF1* zenith = new(angle[i]) TF1("zenith",zenithalDisributionFunction, fZenithMin, fZenithMax);
503
504         // Fix the momentum dependent coefficients
505         zenith->SetParName(0, "a(p)");
506         zenith->SetParameter(0, fAp->At(i));
507
508       }
509
510     } 
511
512   }
513
514 }
515
516 //____________________________________________________________________________
517 const Float_t AliGenCRT::GetZenithAngle(Float_t mom)
518 {
519
520   Float_t zenith = 0.;
521   // Check if you need to generate a constant zenith angle.
522   if ( !fZenithDist ) {
523     // Check if you have defined an array of momentum functions
524     if ( fPDist ) {
525       Int_t pIndex = TMath::Abs(TMath::Nint(mom));
526       TF1* zenithAngle = (TF1*)fPDist->UncheckedAt(pIndex);
527       zenith = zenithAngle->GetRandom();
528       return zenith;
529     } else {
530
531       if ( fCRMode != kMuonFlux ) {
532         // If you aren't generating muons obeying any ditribution
533         // only generate a flat zenith angle, acording the input settings
534         Float_t random[2];
535         Rndm(random, 2);
536         zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
537
538       } else {
539         // Even if you are generating muons acording some distribution,
540         // but you don't want to ...
541         zenith = fZenithMin;
542
543       }
544
545     }
546   } else {
547     zenith = fZenithDist->GetRandom();
548   }
549      
550   return zenith;
551   
552 }
553
554 //____________________________________________________________________________