]>
Commit | Line | Data |
---|---|---|
550bfb3a | 1 | |
2 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
3 | ||
4 | #include "TH1.h" | |
5 | #include "TFile.h" | |
6 | #include "TCanvas.h" | |
7 | #include "TProfile.h" | |
8 | #include "TClonesArray.h" | |
9 | ||
10 | #include "AliModule.h" | |
11 | #include "AliFRAMEv1.h" | |
12 | ||
13 | #include "AliTRDv1.h" | |
14 | #include "AliTRDhit.h" | |
15 | #include "AliTRDsim.h" | |
16 | #include "AliTRDsimple.h" | |
17 | #include "AliTRDsimpleGen.h" | |
18 | #include "AliTRDsimpleMC.h" | |
19 | #include "AliTRDgeometry.h" | |
20 | #include "AliTRDdigitizer.h" | |
21 | #include "AliTRDdigitsManager.h" | |
22 | #include "AliTRDdataArrayI.h" | |
23 | ||
24 | #endif | |
25 | ||
26 | void simple() | |
27 | { | |
28 | ||
29 | /////////////////////////////////////////////////////////////////////////////// | |
30 | // | |
31 | // Macro to run the simplified simulator | |
32 | // | |
33 | /////////////////////////////////////////////////////////////////////////////// | |
34 | ||
35 | //_____________________________________________________________________________ | |
36 | // | |
37 | // Initialization part | |
38 | //_____________________________________________________________________________ | |
39 | // | |
40 | ||
41 | // The simple simulator | |
42 | AliTRDsimple *simple = new AliTRDsimple(); | |
43 | ||
44 | // Initialize a dummy frame so that the TRD works | |
45 | new AliFRAMEv1("FRAME","Space Frame"); | |
46 | ||
47 | // Initialize the TRD detector | |
48 | AliTRDv1 *trd = new AliTRDv1("TRD","TRD slow simulator"); | |
49 | ||
50 | // Needed for some material properties | |
51 | trd->CreateMaterials(); | |
52 | ||
53 | // Use PAI model (most propable energy loss) for the relativistic rise | |
54 | //trd->SetPAI(1); | |
55 | // Use GEANT (mean energy loss) for the relativistic rise | |
56 | trd->SetPAI(0); | |
57 | ||
58 | // Use the Ermilova spectrum for the delta electrons | |
59 | trd->SetErmilova(1); | |
60 | // Use the GEANT spectrum for the delta electrons | |
61 | //trd->SetErmilova(0); | |
62 | ||
63 | // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 85% Xe + 15% CO2) | |
64 | trd->SetGasMix(1); | |
65 | ||
66 | // Get the pointer to the geometry object | |
67 | AliTRDgeometry *geometry = trd->GetGeometry(); | |
68 | ||
69 | // The number of timebins | |
70 | geometry->SetNTimeBin(30); | |
71 | ||
72 | // The additional timebins before and after the drift region | |
73 | geometry->SetExpandTimeBin(5,15); | |
74 | ||
75 | // Switch on TR production | |
76 | AliTRDsim *trdsim = trd->CreateTR(); | |
77 | trdsim->SetNFoils(100); | |
78 | trdsim->SetFoilThick(0.0013); | |
79 | trdsim->SetGapThick(0.0500); | |
80 | ||
81 | // Initialize the TRD object | |
82 | trd->Init(); | |
83 | ||
84 | // The digitizer | |
85 | AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer","Digitizer class"); | |
86 | digitizer->SetSimple(); | |
87 | digitizer->SetGasGain(3200.); // GEANT + new Ermilova spectrum | |
88 | digitizer->SetChipGain(6.1); | |
89 | digitizer->SetNoise(1000.); | |
90 | digitizer->SetADCinRange(1000.); | |
91 | digitizer->SetADCoutRange(1023.); | |
92 | digitizer->SetADCthreshold(0); | |
93 | digitizer->InitDetector(); | |
94 | //digitizer->SetTimeResponse(0); | |
95 | //digitizer->SetVerbose(1); | |
96 | ||
97 | // The event generator | |
98 | AliTRDsimpleGen *generator = simple->GetGenerator(); | |
99 | generator->SetMomentum(3.0,3.0); | |
100 | generator->SetPdg(11); // Electrons | |
101 | //generator->SetPdg(211); // Pions | |
102 | ||
103 | //_____________________________________________________________________________ | |
104 | // | |
105 | // Histograms | |
106 | //_____________________________________________________________________________ | |
107 | // | |
108 | ||
109 | Int_t timeMax = geometry->GetTimeTotal(); | |
110 | Int_t adcRange = ((Int_t) digitizer->GetADCoutRange()); | |
111 | ||
112 | TH1F *hQ = new TH1F("hQ" ,"Charge per hit (all)" ,100,0.0,1000.0); | |
113 | TH1F *hQdedx = new TH1F("hQdedx","Charge per hit (dedx)",100,0.0,1000.0); | |
114 | TH1F *hQtr = new TH1F("hQtr ","Charge per hit (tr)" ,100,0.0,1000.0); | |
115 | ||
116 | TProfile *hQX = new TProfile("hQX" ,"Charge per hit vs X (all)" ,35,0.0,3.5); | |
117 | TProfile *hQXdedx = new TProfile("hQXdedx","Charge per hit vs X (dedx)",35,0.0,3.5); | |
118 | TProfile *hQXtr = new TProfile("hQXtr ","Charge per hit vs X (tr)" ,35,0.0,3.5); | |
119 | ||
120 | TH1F *hNstep = new TH1F("hNstep","Number of steps / cm" , 151,-0.5, 150.5); | |
121 | TH1F *hNel = new TH1F("hNel" ,"Number of electrons / cm",1001,-0.5,1000.5); | |
122 | ||
123 | TH1F *hAmp = new TH1F("hAmp","Amplitude of the digits" | |
124 | ,adcRange+1,-0.5,((Float_t) adcRange)+0.5); | |
125 | TProfile *hAmpTime = new TProfile("hAmpTime","Amplitude vs timebin" | |
126 | ,timeMax,-0.5,((Float_t) timeMax)-0.5); | |
127 | ||
128 | //_____________________________________________________________________________ | |
129 | // | |
130 | // Event loop | |
131 | //_____________________________________________________________________________ | |
132 | // | |
133 | ||
134 | // Number of events (tracks) | |
135 | Int_t nEvent = 10000; | |
136 | ||
137 | Float_t x0 = geometry->GetTime0(0) - AliTRDgeometry::DrThick(); | |
138 | ||
139 | TClonesArray *hitsArray = trd->Hits(); | |
140 | ||
141 | for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) { | |
142 | ||
143 | if (!(iEvent % 100) && (iEvent)) { | |
144 | printf("Event no. %d\n",iEvent); | |
145 | } | |
146 | ||
147 | // Generate the hits for one track | |
148 | simple->ProcessEvent(iEvent); | |
149 | ||
150 | // Analyze the hits | |
151 | Float_t nElectrons = 0.0; | |
152 | for (Int_t iHit = 0; iHit < hitsArray->GetEntries(); iHit++) { | |
153 | ||
154 | AliTRDhit *hit = (AliTRDhit *) hitsArray->UncheckedAt(iHit); | |
155 | Int_t charge = TMath::Abs(hit->GetCharge()); | |
156 | Float_t x = hit->X() - x0; | |
157 | hQ->Fill(charge); | |
158 | hQX->Fill(x,charge); | |
159 | if (hit->FromDrift() || | |
160 | hit->FromAmplification()) { | |
161 | hQdedx->Fill(charge); | |
162 | hQXdedx->Fill(x,charge); | |
163 | nElectrons += charge; | |
164 | } | |
165 | if (hit->FromTRphoton()) { | |
166 | hQtr->Fill(charge); | |
167 | hQXtr->Fill(x,charge); | |
168 | } | |
169 | ||
170 | } | |
171 | ||
172 | hNstep->Fill(((AliTRDsimpleMC *) gMC)->GetNStep() / 3.5); | |
173 | hNel->Fill(nElectrons / 3.5); | |
174 | ||
175 | // Digitize the hits | |
176 | digitizer->MakeDigits(); | |
177 | ||
178 | // Analyze the digits | |
179 | Int_t row = 8; | |
180 | Int_t col = 48; | |
181 | AliTRDdigitsManager *digitsManager = digitizer->Digits(); | |
182 | AliTRDdataArrayI *digitsArray = digitsManager->GetDigits(12); | |
183 | for (Int_t time = 0; time < timeMax; time++) { | |
184 | ||
185 | Int_t amp = digitsArray->GetDataUnchecked(row,col,time); | |
186 | if (amp > 0) { | |
187 | hAmp->Fill((Double_t) amp); | |
188 | hAmpTime->Fill((Double_t) time, (Double_t) amp); | |
189 | } | |
190 | ||
191 | } | |
192 | ||
193 | } | |
194 | ||
195 | //_____________________________________________________________________________ | |
196 | // | |
197 | // Plot the spectra | |
198 | //_____________________________________________________________________________ | |
199 | // | |
200 | ||
201 | TCanvas *cHit = new TCanvas("cHit","Hits",50,50,800,600); | |
202 | cHit->Divide(3,2); | |
203 | ||
204 | cHit->cd(1); | |
205 | gPad->SetLogy(); | |
206 | hQ->SetLineColor(4); | |
207 | hQ->SetXTitle("Q (a.u.)"); | |
208 | hQ->SetYTitle("Entries"); | |
209 | hQ->Draw(); | |
210 | hQdedx->SetLineColor(3); | |
211 | hQdedx->Draw("SAME"); | |
212 | hQtr->SetLineColor(2); | |
213 | hQtr->Draw("SAME"); | |
214 | ||
215 | cHit->cd(2); | |
216 | gPad->SetLogy(); | |
217 | hQX->SetLineColor(4); | |
218 | hQX->SetXTitle("x (cm)"); | |
219 | hQX->SetYTitle("<Q> (a.u.)"); | |
220 | hQX->Draw("HIST"); | |
221 | hQXdedx->SetLineColor(3); | |
222 | hQXdedx->Draw("SAMEHIST"); | |
223 | ||
224 | cHit->cd(3); | |
225 | gPad->SetLogy(); | |
226 | hQXtr->SetLineColor(2); | |
227 | hQXtr->SetXTitle("x (cm)"); | |
228 | hQXtr->SetYTitle("<Q> (a.u.)"); | |
229 | hQXtr->Draw("HIST"); | |
230 | ||
231 | cHit->cd(4); | |
232 | hNstep->SetLineColor(4); | |
233 | hNstep->SetXTitle("N_{step} / cm"); | |
234 | hNstep->SetYTitle("Entries"); | |
235 | hNstep->Draw(); | |
236 | ||
237 | cHit->cd(5); | |
238 | hNel->SetLineColor(4); | |
239 | hNel->SetXTitle("N_{electron} / cm"); | |
240 | hNel->SetYTitle("Entries"); | |
241 | hNel->Draw(); | |
242 | ||
243 | TCanvas *cDigit = new TCanvas("cDigit","Digits",50,50,600,400); | |
244 | cDigit->Divide(2,1); | |
245 | ||
246 | cDigit->cd(1); | |
247 | gPad->SetLogy(); | |
248 | hAmp->SetLineColor(2); | |
249 | hAmp->SetXTitle("ADC channels"); | |
250 | hAmp->SetYTitle("Entries"); | |
251 | hAmp->Draw(); | |
252 | ||
253 | cDigit->cd(2); | |
254 | hAmpTime->SetLineColor(2); | |
255 | hAmpTime->SetXTitle("Timebin"); | |
256 | hAmpTime->SetYTitle("<ADC channels>"); | |
257 | hAmpTime->Draw("HIST"); | |
258 | ||
259 | //_____________________________________________________________________________ | |
260 | // | |
261 | // Save the histograms | |
262 | //_____________________________________________________________________________ | |
263 | // | |
264 | ||
265 | TFile *fileOut = new TFile("simple.root","RECREATE"); | |
266 | ||
267 | hQ->Write(); | |
268 | hQdedx->Write(); | |
269 | hQtr->Write(); | |
270 | ||
271 | hQX->Write(); | |
272 | hQXdedx->Write(); | |
273 | hQXtr->Write(); | |
274 | ||
275 | hNstep->Write(); | |
276 | hNel->Write(); | |
277 | ||
278 | hAmp->Write(); | |
279 | hAmpTime->Write(); | |
280 | ||
281 | fileOut->Close(); | |
282 | ||
283 | } |