cbd85fb45234524cfbac581fded5a960cedbcb34
[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 /* $Id$ */
17
18 /////////////////////////////////////////////////////////////////////////////
19 //
20 //  Contain parametrizations to generate atmospheric muons, and also
21 //  to generate single muons and muon bundles at surface level.
22 //
23 //Begin_Html
24 /*
25 <img src="picts/AliGenCRTClass.gif">
26 </pre>
27 <br clear=left>
28 <font size=+2 color=red>
29 <p>The responsible person for this module is
30 <a href="mailto:Enrique.Gamez.Flores@cern.ch">Enrique Gamez</a>.
31 </font>
32 <pre>
33 */
34 //End_Html
35 //
36 /////////////////////////////////////////////////////////////////////////////
37
38 #include "AliGenCRT.h"
39
40 #include <TMCProcess.h>
41 #include <TPDGCode.h>
42 #include <TClonesArray.h>
43 #include <TF1.h>
44
45 #include "AliRun.h"
46 #include "AliConst.h"
47
48 ClassImp(AliGenCRT)
49
50 //_____________________________________________________________________________
51 AliGenCRT::AliGenCRT()
52   : AliGenerator(),
53     fIpart(0),
54     fCRMode(kSingleMuons),
55     fCRModeName(0),
56     fXwidth(0),
57     fNx(1),
58     fZwidth(0),
59     fNz(1),
60     fMuonGrid(kFALSE),
61     fZenithMin(0),
62     fZenithMax(0),
63     fAzimuthMin(0),
64     fAzimuthMax(0),
65     fPRange(0),
66     fPResolution(1),
67     fAp(0),
68     fMomentumDist(0),
69     fUnfoldedMomentumDist(0),
70     fZenithDist(0),
71     fPDist(0)
72 {
73   //
74   // Default ctor.
75   //
76 }
77
78 //_____________________________________________________________________________
79 AliGenCRT::AliGenCRT(Int_t npart) 
80   : AliGenerator(npart),
81     fIpart(kMuonMinus),
82     fCRMode(kSingleMuons),
83     fCRModeName(0),
84     fXwidth(0),
85     fNx(1),
86     fZwidth(0),
87     fNz(1),
88     fMuonGrid(kFALSE),
89     fZenithMin(0),
90     fZenithMax(0),
91     fAzimuthMin(0),
92     fAzimuthMax(0),
93     fPRange(0),
94     fPResolution(1),
95     fAp(0),
96     fMomentumDist(0),
97     fUnfoldedMomentumDist(0),
98     fZenithDist(0),
99     fPDist(0)
100 {
101   //
102   // Standard ctor.
103   //
104   fName = "CRT";
105   fTitle = "Cosmic Muons generator";
106
107   // Set the origin above the vertex, on the surface.
108   fOrigin[0] = 0.;
109   fOrigin[1] = AliCRTConstants::Instance()->Depth(); // At the surface by default.
110   fOrigin[2] = 0.;
111 }
112
113 //_____________________________________________________________________________
114 AliGenCRT::AliGenCRT(const AliGenCRT& gen)
115   : AliGenerator(gen)
116 {
117   //
118   // Copy constructor
119   //
120   gen.Copy(*this);
121 }
122
123 //_____________________________________________________________________________
124 AliGenCRT& AliGenCRT::operator=(const AliGenCRT& gen)
125 {
126   //
127   // Asingment operator
128   //
129   gen.Copy(*this);
130   return *this;
131 }
132
133 //_____________________________________________________________________________
134 AliGenCRT::~AliGenCRT()
135 {
136   //
137   // Default dtor.
138   //
139   if ( fPDist ) {fPDist->Delete(); delete fPDist; fPDist = 0;}
140   if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
141   if ( fMomentumDist ) delete fMomentumDist;
142   if ( fAp )           delete fAp;
143   if ( fCRModeName )   delete fCRModeName;
144 }
145
146 //_____________________________________________________________________________
147 void AliGenCRT::Generate()
148 {
149   //
150   // Generate on one trigger
151   // Call the respective method inside the loop for the number
152   // of tracks per trigger.
153
154   for (Int_t i = 0; i < fNpart; i++ ) {
155
156     if ( fCRMode == kMuonBundle ) {
157       this->GenerateOneMuonBundle();
158
159     } else if ( fCRMode == kSingleMuons ) {
160       this->GenerateOneSingleMuon(kTRUE);
161
162     } else {
163       // Generate only single muons but following the parametrizations
164       // for atmospheric muons.
165       this->GenerateOneSingleMuon();
166
167     }
168
169   }
170 }
171
172 //_____________________________________________________________________________
173 void AliGenCRT::Init()
174 {
175   //
176   // Initialize some internal methods.
177   //
178
179   // Determine some specific data members.
180   fPRange = TMath::Abs(fPMax-fPMin);
181
182   if ( fCRMode == kSingleMuons ) {
183     fCRModeName = new TString("Single Muons");
184     // Initialisation, check consistency of selected ranges
185     if(TestBit(kPtRange)&&TestBit(kMomentumRange)) 
186       Fatal("Init","You should not set the momentum range and the pt range!");
187     
188     if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange))) 
189       Fatal("Init","You should set either the momentum or the pt range!");
190     
191   } else if ( fCRMode == kMuonBundle ) {
192     fCRModeName = new TString("Muon Bundles");
193
194   } else if ( fCRMode == kMuonFlux ) {
195     fCRModeName = new TString("Muon Fluxes");
196     // Initialize the ditribution functions.
197     this->InitMomentumGeneration();
198     this->InitZenithalAngleGeneration();
199     
200   } else {
201     Fatal("Generate", "Generation Mode unknown!\n");
202
203   }
204
205 }
206
207 //____________________________________________________________________________
208 void AliGenCRT::GenerateOneSingleMuon(Bool_t withFlatMomentum)
209 {
210   //
211   // Generate One Single Muon
212   // This method will generate only one muon.
213   // The momentum will be randomly flat distributed if
214   // the paremeter "withFlatMomentum" is set to kTRUE,
215   // otherwise the momentum will generate acordingly the parametrization
216   // given by 
217   // and adpted from Bruno Alessandro's implementation with the
218   // CERNLIB to AliRoot.
219   // The "withFlatMomentum" parameter also will be used to generate
220   // the muons with a flat Zenithal angle distribution.
221   // Do the smearing here, so that means per track.
222
223   Float_t polar[3]= {0,0,0}; // Polarization parameters
224   Float_t origin[3];
225   Int_t nt;
226   Float_t p[3];
227   Float_t pmom, pt;
228   Float_t random[6];
229
230   // Take the azimuth random.
231   Rndm(random, 2);
232   Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
233   Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
234
235   if ( withFlatMomentum ) {
236     Rndm(random,3);
237     if(TestBit(kMomentumRange)) {
238       pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
239       pt = pmom*TMath::Sin(zenith*kDegrad);
240     } else {
241       pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
242       pmom = pt/TMath::Sin(zenith*kDegrad);
243     }
244
245   } else {
246     if ( fMomentumDist ) {
247       pmom = -this->GetMomentum(); // Always downwards.
248     } else {
249       pmom = -fPMin;
250     }
251     zenith = this->GetZenithAngle(pmom);  // In degrees
252     pt = pmom*TMath::Sin(zenith*kDegrad);
253   }
254
255   p[0] = pt*TMath::Sin(azimuth*kDegrad);
256   p[1] = pmom*TMath::Cos(zenith*kDegrad);
257   p[2] = pt*TMath::Cos(azimuth*kDegrad);
258
259   // Finaly the origin, with the smearing
260   Rndm(random,6);
261   origin[0] = AliCRTConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
262     TMath::Sin(azimuth*kDegrad);
263   // + fOsigma[0]* TMath::Cos(2*random[0]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[1]));
264
265   origin[1] = AliCRTConstants::Instance()->Depth();
266
267   origin[2] = AliCRTConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
268     TMath::Cos(azimuth*kDegrad);
269   // + fOsigma[2]* TMath::Cos(2*random[2]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[3]));;
270
271   // Put the track on the stack.
272   PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
273
274 }
275
276 //____________________________________________________________________________
277 void AliGenCRT::GenerateOneMuonBundle()
278 {
279   //
280   // Generate One Muon Bundle method
281   // This method will generate a bunch of muons following the
282   // procedure of the AliGenScan class.
283   // These muons will be generated with flat momentum.
284
285   Float_t polar[3]= {0,0,0}; // Polarization parameters
286   Float_t origin[3];
287   Float_t p[3];
288   Int_t nt;
289   Float_t pmom;
290   Float_t random[6];
291
292   Rndm(random, 3);
293   //Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
294   //Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[2];
295   Float_t zenith = 10;
296   Float_t azimuth = 30;
297
298   // Generate the kinematics a la AliGenScan (Andreas Morchs)
299   Float_t dx, dz;
300   if ( fNx > 0 ) {
301     dx = fXwidth/fNx;
302   } else {
303     dx = 1e10;
304     //dx = 100.;
305   }
306
307   if ( fNz > 0 ) {
308     dz = fZwidth/fNz;
309   } else {
310     dz = 1e10;
311     //dz = 100.;
312   }
313
314   //origin[0] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
315   origin[0] = 0.;
316     //TMath::Sin(azimuth*kDegrad);
317   //origin[1] = AliCRTConstants::fgDepth;
318   origin[1] = 900;
319   //origin[2] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
320   origin[2] = 0.;
321     //TMath::Cos(azimuth*kDegrad);
322
323   for (Int_t ix = 0; ix < fNx; ix++ ) {
324     for (Int_t iz = 0; iz < fNz; iz++ ) {
325       Rndm(random,6);
326       origin[0]+=ix*dx+2*(random[1]-0.5)*fOsigma[0];
327       origin[2]+=iz*dz+2*(random[2]-0.5)*fOsigma[2];
328       if ( random[4] < 0.5 ) {
329         origin[0] = -1*origin[0];
330       }
331       if ( random[5] < 0.5 ) {
332         origin[2] = -1*origin[2];
333       }
334
335       pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downwards
336       p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
337       p[1] = TMath::Cos(zenith*kDegrad)*pmom;
338       p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
339
340       PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
341     }
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       this->InitApWeightFactors();
487       
488       // Define the standard function.
489       const 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) const
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
522       Float_t tmpzenith = TMath::Cos(zenithAngle->GetRandom()*kDegrad);
523       // Correct the value
524       zenith = kRaddeg*TMath::ACos(1 - (tmpzenith - 1)/fAp->At(pIndex));
525       return zenith;
526     } else {
527
528       if ( fCRMode != kMuonFlux ) {
529         // If you aren't generating muons obeying any ditribution
530         // only generate a flat zenith angle, acording the input settings
531         Float_t random[2];
532         Rndm(random, 2);
533         zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
534
535       } else {
536         // Even if you are generating muons acording some distribution,
537         // but you don't want to ...
538         zenith = fZenithMin;
539
540       }
541
542     }
543   } else {
544     zenith = fZenithDist->GetRandom();
545   }
546
547   return zenith;
548 }
549
550 //_____________________________________________________________________________
551 const Float_t AliGenCRT::GetMomentum() const
552 {
553   //
554   //
555   //
556   return fMomentumDist->GetRandom();
557 }