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