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