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