Update of simple generator
[u/mrichter/AliRoot.git] / TRD / AliTRDrunSimple.C
CommitLineData
fa148e6c 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
26void AliTRDrunSimple()
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 simple->Init();
44
45 // Initialize a dummy frame so that the TRD works
46 new AliFRAMEv1("FRAME","Space Frame");
47
48 // Initialize the TRD detector
49 AliTRDv1 *trd = new AliTRDv1("TRD","TRD slow simulator");
ee7e3ede 50 trd->SetHitTypeStandard();
fa148e6c 51
52 // Needed for some material properties
53 trd->CreateMaterials();
54
55 // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 85% Xe + 15% CO2)
56 trd->SetGasMix(1);
57
ee7e3ede 58 // Define the parameter object
59 // If no external parameter object is defined,
60 // default parameter will be used
61 AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
62 ,"TRD parameter class");
63
64 parameter->SetADCthreshold(0);
65 parameter->SetNTimeBin(30); // The number of timebins
66 parameter->SetExpandTimeBin(5,15); // The additional timebins
fa148e6c 67
68 // Switch on TR production
69 AliTRDsim *trdsim = trd->CreateTR();
70 trdsim->SetNFoils(100);
71 trdsim->SetFoilThick(0.0013);
72 trdsim->SetGapThick(0.0500);
73
74 // Initialize the TRD object
75 trd->Init();
76
77 // The digitizer
78 AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer","Digitizer class");
ee7e3ede 79
80 digitizer->SetParameter(parameter);
fa148e6c 81 digitizer->SetSimple();
fa148e6c 82 digitizer->InitDetector();
fa148e6c 83
84 // The event generator
85 AliTRDsimpleGen *generator = simple->GetGenerator();
86 generator->SetMomentum(3.0,3.0);
87 generator->SetPdg(11); // Electrons
fa148e6c 88
89 //_____________________________________________________________________________
90 //
91 // Histograms
92 //_____________________________________________________________________________
93 //
94
ee7e3ede 95 Int_t timeMax = parameter->GetTimeTotal();
96 Int_t adcRange = ((Int_t) parameter->GetADCoutRange());
fa148e6c 97
98 TH1F *hQ = new TH1F("hQ" ,"Charge per hit (all)" ,100,0.0,1000.0);
99 TH1F *hQdedx = new TH1F("hQdedx","Charge per hit (dedx)",100,0.0,1000.0);
100 TH1F *hQtr = new TH1F("hQtr ","Charge per hit (tr)" ,100,0.0,1000.0);
101
102 TProfile *hQX = new TProfile("hQX" ,"Charge per hit vs X (all)" ,35,0.0,3.5);
103 TProfile *hQXdedx = new TProfile("hQXdedx","Charge per hit vs X (dedx)",35,0.0,3.5);
104 TProfile *hQXtr = new TProfile("hQXtr ","Charge per hit vs X (tr)" ,35,0.0,3.5);
105
106 TH1F *hNstep = new TH1F("hNstep","Number of steps / cm" , 151,-0.5, 150.5);
107 TH1F *hNel = new TH1F("hNel" ,"Number of electrons / cm",1001,-0.5,1000.5);
108
109 TH1F *hAmp = new TH1F("hAmp","Amplitude of the digits"
110 ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
111 TProfile *hAmpTime = new TProfile("hAmpTime","Amplitude vs timebin"
112 ,timeMax,-0.5,((Float_t) timeMax)-0.5);
113
114 //_____________________________________________________________________________
115 //
116 // Event loop
117 //_____________________________________________________________________________
118 //
119
120 // Number of events (tracks)
121 Int_t nEvent = 10000;
122
ee7e3ede 123 Float_t x0 = parameter->GetTime0(0) - AliTRDgeometry::DrThick();
fa148e6c 124
125 TClonesArray *hitsArray = trd->Hits();
126
127 for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
128
129 if (!(iEvent % 100) && (iEvent)) {
130 printf("Event no. %d\n",iEvent);
131 }
132
133 // Generate the hits for one track
134 simple->ProcessEvent(iEvent);
135
136 // Analyze the hits
137 Float_t nElectrons = 0.0;
138 for (Int_t iHit = 0; iHit < hitsArray->GetEntries(); iHit++) {
139
140 AliTRDhit *hit = (AliTRDhit *) hitsArray->UncheckedAt(iHit);
141 Int_t charge = TMath::Abs(hit->GetCharge());
142 Float_t x = hit->X() - x0;
143 hQ->Fill(charge);
144 hQX->Fill(x,charge);
145 if (hit->FromDrift() ||
146 hit->FromAmplification()) {
147 hQdedx->Fill(charge);
148 hQXdedx->Fill(x,charge);
149 nElectrons += charge;
150 }
151 if (hit->FromTRphoton()) {
152 hQtr->Fill(charge);
153 hQXtr->Fill(x,charge);
154 }
155
156 }
157
158 hNstep->Fill(((AliTRDsimpleMC *) gMC)->GetNStep() / 3.5);
159 hNel->Fill(nElectrons / 3.5);
160
161 // Digitize the hits
162 digitizer->MakeDigits();
163
164 // Analyze the digits
165 Int_t row = 8;
166 Int_t col = 48;
167 AliTRDdigitsManager *digitsManager = digitizer->Digits();
168 AliTRDdataArrayI *digitsArray = digitsManager->GetDigits(12);
169 for (Int_t time = 0; time < timeMax; time++) {
170
171 Int_t amp = digitsArray->GetDataUnchecked(row,col,time);
172 if (amp > 0) {
173 hAmp->Fill((Double_t) amp);
174 hAmpTime->Fill((Double_t) time, (Double_t) amp);
175 }
176
177 }
178
179 }
180
181 //_____________________________________________________________________________
182 //
183 // Plot the spectra
184 //_____________________________________________________________________________
185 //
186
187 TCanvas *cHit = new TCanvas("cHit","Hits",50,50,800,600);
188 cHit->Divide(3,2);
189
190 cHit->cd(1);
191 gPad->SetLogy();
192 hQ->SetLineColor(4);
193 hQ->SetXTitle("Q (a.u.)");
194 hQ->SetYTitle("Entries");
195 hQ->Draw();
196 hQdedx->SetLineColor(3);
197 hQdedx->Draw("SAME");
198 hQtr->SetLineColor(2);
199 hQtr->Draw("SAME");
200
201 cHit->cd(2);
202 gPad->SetLogy();
203 hQX->SetLineColor(4);
204 hQX->SetXTitle("x (cm)");
205 hQX->SetYTitle("<Q> (a.u.)");
206 hQX->Draw("HIST");
207 hQXdedx->SetLineColor(3);
208 hQXdedx->Draw("SAMEHIST");
209
210 cHit->cd(3);
211 gPad->SetLogy();
212 hQXtr->SetLineColor(2);
213 hQXtr->SetXTitle("x (cm)");
214 hQXtr->SetYTitle("<Q> (a.u.)");
215 hQXtr->Draw("HIST");
216
217 cHit->cd(4);
218 hNstep->SetLineColor(4);
219 hNstep->SetXTitle("N_{step} / cm");
220 hNstep->SetYTitle("Entries");
221 hNstep->Draw();
222
223 cHit->cd(5);
224 hNel->SetLineColor(4);
225 hNel->SetXTitle("N_{electron} / cm");
226 hNel->SetYTitle("Entries");
227 hNel->Draw();
228
229 TCanvas *cDigit = new TCanvas("cDigit","Digits",50,50,600,400);
230 cDigit->Divide(2,1);
231
232 cDigit->cd(1);
233 gPad->SetLogy();
234 hAmp->SetLineColor(2);
235 hAmp->SetXTitle("ADC channels");
236 hAmp->SetYTitle("Entries");
237 hAmp->Draw();
238
239 cDigit->cd(2);
240 hAmpTime->SetLineColor(2);
241 hAmpTime->SetXTitle("Timebin");
242 hAmpTime->SetYTitle("<ADC channels>");
243 hAmpTime->Draw("HIST");
244
245 //_____________________________________________________________________________
246 //
247 // Save the histograms
248 //_____________________________________________________________________________
249 //
250
251 TFile *fileOut = new TFile("simple.root","RECREATE");
252
253 hQ->Write();
254 hQdedx->Write();
255 hQtr->Write();
256
257 hQX->Write();
258 hQXdedx->Write();
259 hQXtr->Write();
260
261 hNstep->Write();
262 hNel->Write();
263
264 hAmp->Write();
265 hAmpTime->Write();
266
267 fileOut->Close();
268
269}