StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
frankSimu.C
1 #include <TH2.h>
2 #include <TStyle.h>
3 #include <TCanvas.h>
4 #include <TF1.h>
5 #include <TSystem.h>
6 #include <TLatex.h>
7 #include <iostream>
8 #include <fstream>
9 #include <TGraph.h>
10 #include <TGraphErrors.h>
11 #include <TMath.h>
12 #include <TLine.h>
13 #include <TChain.h>
14 #include <TClonesArray.h>
15 #include <TFile.h>
16 #include <TPostScript.h>
17 #include <TNtuple.h>
18 #include <TString.h>
19 #include <TRandom.h>
20 #include <TRandom2.h>
21 
22 void
23 frankSimu(float rotAngle=0., long nEvents = 1000) {
24 /*
25 gStyle->SetPalette(1);
26 gStyle->SetOptStat(1);
27 gStyle->SetOptFit(11);
28 //gStyle->SetOptTitle(0);
29 gStyle->SetStatBorderSize(1);
30 gStyle->SetStatColor(10);
31 gStyle->SetCanvasColor(10);
32 gStyle->SetPadLeftMargin(0.16);
33 gStyle->SetPadBottomMargin(0.16);
34 gStyle->SetPadTickX(1);
35 gStyle->SetPadTickY(1);
36 gStyle->SetOptTitle(0);
37 gStyle->SetTitleSize(0.048,"xy");
38 gStyle->SetLabelSize(0.04,"xy");
39 gStyle->SetTitleOffset(1.3,"x");
40 gStyle->SetTitleOffset(1.3,"y");
41 */
42 
43 
44 
45 
46  // in order to simulate the energy loss for individual ionizing collisions a table from Bichsel is used.
47  // This is given in figure 9 in NIM A562, 154 (2006)
48  // This table gives the energy loss for a given random number (from 0 to 1 in steps of 10^-4)
49  // Two files, low and high beta gamma:
50  // file Low: beta gamma .31623 1.00000 3.16228 10.00000 31.62278
51  // file High: beta gamma 100.00000 316.22777 1000.00000 3162.27766 10000.00000
52 
53  // here use column 2 for High file
54 
55  double eLoss[10000];
56  ifstream inFile("BichselELossProbHighBG.dat");
57  double cl1, cl2, cl3, cl4, cl5, cl6, cl7;
58  for (int i = 0; i < 10000; i++) {
59  inFile >> cl1 >> cl2 >> cl3 >> cl4 >> cl5 >> cl6>> cl7;
60  eLoss[i] = cl4;
61  }
62 
63  float angle = TMath::Pi() * rotAngle / 180.;
64 
65  double pathLength = 3.2 / TMath::Cos(angle);
66  double pairsPerMM = 4.;
67 
68 
69  TRandom2* tR = new TRandom2();
70  tR->SetSeed(0);
71 
72  TH1F* hMean = new TH1F("hMean", "hMean", 100, -1, 1);
73  TH1F* hNP = new TH1F("hNP", "hNP", 35, -0.5, 34.5);
74  hNP->SetXTitle("Number of Primary Pairs");
75  TH1F* hElectron = new TH1F("hElectron", "hElectron", 250, -0.5, 249.5);
76  TH1F* hEnergy = new TH1F("hEnergy", "hEnergy", 250, 0, 5.);
77  hEnergy->SetXTitle("Energy Loss [keV]");
78  TH1F* hEPerColl = new TH1F("hEPerColl", "hEPerColl", 100, 0, 100);
79  hEPerColl->SetXTitle("Energy Loss per collision [eV]");
80 
81  TF1* fElectronDist = new TF1("fElectronDist", "1/(x*x*x)", 0, 100);
82  float pos[1000];
83  Double_t weight[1000];
84  int totalElectrons = 0;
85  float totalEnergy = 0;
86  double dist = 0;
87 
88  for (int i = 0; i < nEvents; i++) {
89 
90  int np = 0;
91  totalElectrons = 0;
92  totalEnergy = 0;
93  dist = 0;
94  while (1) {
95  // make a random step
96  double stepLength = - TMath::Log(tR->Uniform()) / pairsPerMM;
97  dist += stepLength;
98  if (dist > pathLength) break;
99  np++;
100  pos[np-1] = dist/pathLength -0.5;//tR->Uniform() - 0.5;
101  // additional weight according to secondary energy distribution, according to Bichsel dist
102  int rndBin = ((int) (10000.0 * tR->Uniform()));
103  double eL = eLoss[rndBin];
104  // alternative approach: subtract Ar binding energy , see how much is left over
105  int ns = 1 + ((int) ((eL-15.4)/26.));
106  totalEnergy += eL;
107  hEPerColl->Fill(eL);
108  if (ns < 0) ns = 0;
109  weight[np-1] = ns;
110  totalElectrons += weight[np-1];
111  }
112  float mpos = TMath::Mean(np, pos, weight);
113  mpos = mpos * 3.2 * TMath::Sin(angle);
114  hMean->Fill(mpos);
115  hNP->Fill(np);
116  hElectron->Fill(totalElectrons);
117  totalEnergy = totalEnergy / 1000.;
118  hEnergy->Fill(totalEnergy);
119  }
120 
121  TCanvas* c1 = new TCanvas("SimuMean", "SimuMean", 600, 450);
122  hMean->Draw();
123  TCanvas* c2 = new TCanvas("NumberOfPairs", "NumberOfPairs", 600, 800);
124  c2->Divide(1, 2);
125  c2->cd(1);
126  hNP->Draw();
127  c2->cd(2);
128  hEPerColl->Draw();
129  TCanvas* c3 = new TCanvas("NumberOfElectrons", "NumberOfElectrons", 600, 800);
130  c3->Divide(1, 2);
131  c3->cd(1);
132  hElectron->Draw();
133  c3->cd(2);
134  hEnergy->Draw();
135 
136 }
Generated by   doxygen 1.8.5