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