16 {
17
18 std::cout << "ElectronGenerator: generating " << m_nevents << " electrons" << std::endl;
19 std::cout << "\t\t" << m_upperbg << "\t" << m_lowerbg << std::endl;
20
21 int run = 1;
22 int nhit;
23 double dedxpub;
24 double dedx;
26 double p;
28 double costh;
29 double db;
30 double dz;
31 double chiPi;
32 double eop;
33
34 double dedx_cur, dedx_res;
35
36
37 TTree* gentrack = new TTree("track","Fake track sample");
38
39 gentrack->Branch("run",&run,"run/I");
40 gentrack->Branch("numGoodHits",&nhit,"numGoodHits/I");
41 gentrack->Branch("dedx",&dedxpub,"dedx/D");
42 gentrack->Branch("dedxcur",&dedx_cur,"dedxcur/D");
43 gentrack->Branch("dedxres",&dedx_res,"dedxres/D");
44 gentrack->Branch("dedxsat",&dedx,"dedxsat/D");
45 gentrack->Branch(
"q",&
q,
"q/D");
46 gentrack->Branch("pF",&p,"pF/D");
47 gentrack->Branch(
"bg",&
bg,
"bg/D");
48 gentrack->Branch("costh",&costh,"costh/D");
49 gentrack->Branch("db",&db,"db/D");
50 gentrack->Branch("dz",&dz,"dz/D");
51 gentrack->Branch("chiPi",&chiPi,"chiPi/D");
52 gentrack->Branch("eopst",&eop,"eopst/D");
53
54
55 TTree* genhit = new TTree("hit","Fake hit sample");
56
57 const int kMaxHits = 100;
58 double doca[kMaxHits];
59 double enta[kMaxHits];
60 double dedxhit[kMaxHits];
61
62 genhit->Branch("numGoodHits",&nhit,"numGoodHits/I");
63 genhit->Branch("doca",doca,"doca[numGoodHits]/D");
64 genhit->Branch("enta",enta,"enta[numGoodHits]/D");
65 genhit->Branch("dedxhit",dedxhit,"dedxhit[numGoodHits]/D");
66 genhit->Branch("dedx",&dedxpub,"dedx/D");
67 genhit->Branch("dedxsat",&dedx,"dedxsat/D");
68 genhit->Branch("pF",&p,"pF/D");
69 genhit->Branch("costh",&costh,"costh/D");
70 genhit->Branch("db",&db,"db/D");
71 genhit->Branch("dz",&dz,"dz/D");
72 genhit->Branch("chiPi",&chiPi,"chiPi/D");
73
74 TH2F* dedx_costh = new TH2F("dedx_costh","dE/dx vs. cos(#theta);cos(#theta);dE/dx",
75 100,-1.0,1.0,100,0.0,10.0);
76 TH2F* dedx_pq = new TH2F("dedx_pq","dE/dx vs. pq;pq [GeV];dE/dx",
77 100,-1.0*m_upperbg,m_upperbg,100,0.8,1.2);
78
79 TRandom *r0 = new TRandom();
80
81 r0->SetSeed(0);
82
83 double mass = Widget::melectron;
84
85 TF1 empcor("empcor","1-0.15*TMath::Exp(-5*x)",0.2,2);
86 TF1 docashape("docashape","-78.3*x^2+75.5",-1,1);
87 TF1 enta1("enta1","1500*TMath::Landau(x,0.125,0.05)",0,3);
88 TF1 enta2("enta2","2000*TMath::Landau(TMath::Abs(x),0.125,0.05)",-3,0);
89 TF1 entashape("entashape","enta1+enta2",-3,3);
90 TF1 satshape("satshape","1.0-2.0*TMath::Landau(TMath::Abs(x),0.01,0.05)",-1,1);
91
92
93 for( int i = 0; i < m_nevents; ++i ){
94 if( i % 100 == 0 ){
95 std::cout << i << " events" << std::endl;
96 run++;
97 }
98
99 double q = (r0->Rndm(i) < 0.5) ? -1.0 : 1.0;
100
101 p = r0->Rndm(i)*(m_upperbg-m_lowerbg)+m_lowerbg;
102 while( p < m_lowerbg || p > m_upperbg )
103 p = r0->Rndm(i)*(m_upperbg-m_lowerbg)+m_lowerbg;
105
106 costh = 2*(r0->Rndm(i))-1;
107 while( costh < -1.0 || costh > 1.0 )
108 costh = 2*(r0->Rndm(i))-1;
109
110 dedx = (
q < 0) ? (2 - empcor.Eval(p)) : empcor.Eval(p);
111 dedx = dedx*satshape.Eval(costh);
112 dedx_cur = 1.0;
113 dedx_res = 1.05-empcor.Eval(p);
114
115 nhit = r0->Gaus(25);
116 while( nhit > 100 || nhit < 0 )
117 nhit = r0->Gaus(25);
118 eop = 0.5;
119 chiPi = r0->Gaus(1.0);
120
121 for( int j = 0; j < nhit; ++j ){
122 enta[j] = entashape.GetRandom();
123 doca[j] = docashape.GetRandom();
124 dedxhit[j] = r0->Landau(dedx,dedx_res);
125 }
126
127 dedxpub = dedx;
128
129 dedx_costh->Fill(costh,dedx);
130 dedx_pq->Fill(p*
q,dedx);
131
132 gentrack->Fill();
133 genhit->Fill();
134 }
135 delete r0;
136
137 genfile->cd();
138 gentrack->Write();
139 genhit->Write();
140 dedx_costh->Write();
141 dedx_pq->Write();
142 genfile->Close();
143}
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
1.9.6