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