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