Trigger board object base class
[u/mrichter/AliRoot.git] / CRT / AliGenCRT.cxx
CommitLineData
71352da6 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
88cb7938 16/* $Id$ */
71352da6 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
fddb5247 38#include "AliGenCRT.h"
39
116cbefd 40#include <TMCProcess.h>
41#include <TPDGCode.h>
bb5ce8b0 42#include <TClonesArray.h>
43#include <TF1.h>
81d50d27 44#include <TH1F.h>
71352da6 45
bb5ce8b0 46#include "AliRun.h"
116cbefd 47#include "AliConst.h"
bb5ce8b0 48
71352da6 49ClassImp(AliGenCRT)
50
51//_____________________________________________________________________________
bb5ce8b0 52AliGenCRT::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)
71352da6 73{
74 //
75 // Default ctor.
76 //
71352da6 77}
78
79//_____________________________________________________________________________
80AliGenCRT::AliGenCRT(Int_t npart)
bb5ce8b0 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)
71352da6 101{
102 //
103 // Standard ctor.
104 //
105 fName = "CRT";
106 fTitle = "Cosmic Muons generator";
107
71352da6 108 // Set the origin above the vertex, on the surface.
109 fOrigin[0] = 0.;
fddb5247 110 fOrigin[1] = AliCRTConstants::Instance()->Depth(); // At the surface by default.
71352da6 111 fOrigin[2] = 0.;
71352da6 112}
113
114//_____________________________________________________________________________
bb5ce8b0 115AliGenCRT::AliGenCRT(const AliGenCRT& gen)
116 : AliGenerator(gen)
71352da6 117{
118 //
bb5ce8b0 119 // Copy constructor
71352da6 120 //
121 gen.Copy(*this);
122}
123
124//_____________________________________________________________________________
bb5ce8b0 125AliGenCRT& AliGenCRT::operator=(const AliGenCRT& gen)
71352da6 126{
127 //
bb5ce8b0 128 // Asingment operator
71352da6 129 //
130 gen.Copy(*this);
131 return *this;
132}
133
134//_____________________________________________________________________________
135AliGenCRT::~AliGenCRT()
136{
137 //
138 // Default dtor.
139 //
bb5ce8b0 140 if ( fPDist ) {fPDist->Delete(); delete fPDist; fPDist = 0;}
71352da6 141 if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
bb5ce8b0 142 if ( fMomentumDist ) delete fMomentumDist;
143 if ( fAp ) delete fAp;
144 if ( fCRModeName ) delete fCRModeName;
71352da6 145}
146
147//_____________________________________________________________________________
148void AliGenCRT::Generate()
149{
150 //
151 // Generate on one trigger
bb5ce8b0 152 // Call the respective method inside the loop for the number
153 // of tracks per trigger.
71352da6 154
bb5ce8b0 155 for (Int_t i = 0; i < fNpart; i++ ) {
71352da6 156
bb5ce8b0 157 if ( fCRMode == kMuonBundle ) {
158 this->GenerateOneMuonBundle();
71352da6 159
bb5ce8b0 160 } else if ( fCRMode == kSingleMuons ) {
161 this->GenerateOneSingleMuon(kTRUE);
71352da6 162
71352da6 163 } else {
81d50d27 164 // Generate only single muons following the parametrizations
bb5ce8b0 165 // for atmospheric muons.
166 this->GenerateOneSingleMuon();
71352da6 167
71352da6 168 }
169
71352da6 170 }
71352da6 171}
172
173//_____________________________________________________________________________
174void 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//____________________________________________________________________________
bb5ce8b0 209void 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);
fddb5247 262 origin[0] = AliCRTConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
a86d3e7e 263 TMath::Sin(azimuth*kDegrad)
81d50d27 264 + fOsigma[0]* TMath::Cos(2*random[0]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[1]));
bb5ce8b0 265
fddb5247 266 origin[1] = AliCRTConstants::Instance()->Depth();
bb5ce8b0 267
fddb5247 268 origin[2] = AliCRTConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
a86d3e7e 269 TMath::Cos(azimuth*kDegrad)
270 + fOsigma[2]* TMath::Cos(2*random[2]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[3]));
bb5ce8b0 271
272 // Put the track on the stack.
273 PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
274
275}
276
277//____________________________________________________________________________
278void 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);
81d50d27 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;
bb5ce8b0 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
81d50d27 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.;
bb5ce8b0 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//____________________________________________________________________________
71352da6 349void 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//____________________________________________________________________________
369void 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//___________________________________________________________________________
405void 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//____________________________________________________________________________
473void 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)
fddb5247 487 this->InitApWeightFactors();
71352da6 488
489 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
81d50d27 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;
71352da6 496 for ( Int_t i = 0; i < pEnd; i++ ) {
497 // Fill the distribution
81d50d27 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 }
71352da6 510
511 }
512
513 }
514
515 }
516
517}
518
519//____________________________________________________________________________
17323043 520Float_t AliGenCRT::GetZenithAngle(Float_t mom) const
71352da6 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));
81d50d27 529 TH1F* cosZenithAngle = (TH1F*)fPDist->UncheckedAt(pIndex);
530 Float_t tmpzenith = TMath::ACos(cosZenithAngle->GetRandom());
bb5ce8b0 531 // Correct the value
81d50d27 532 zenith = kRaddeg*tmpzenith;
71352da6 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 }
bb5ce8b0 554
71352da6 555 return zenith;
71352da6 556}
557
bb5ce8b0 558//_____________________________________________________________________________
17323043 559Float_t AliGenCRT::GetMomentum() const
bb5ce8b0 560{
561 //
562 //
563 //
564 return fMomentumDist->GetRandom();
565}