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