]>
Commit | Line | Data |
---|---|---|
0ac8e4b6 | 1 | |
2 | Int_t AliTRDtest() | |
3 | { | |
4 | // | |
5 | // Test macro for the TRD code | |
6 | // | |
7 | ||
8 | Int_t rc = 0; | |
9 | ||
10 | AliSimulation sim; | |
11 | sim.SetConfigFile("$(ALICE_ROOT)/TRD/Macros/AliTRDconfig.C"); | |
12 | sim.SetLoadAlignFromCDB(0); | |
13 | sim.Run(); | |
14 | ||
15 | // Analyze the TRD hits | |
16 | if (rc = AliTRDanalyzeHits()) return rc; | |
17 | ||
18 | // Analyze the digits | |
19 | //if (rc = AliTRDanalyzeDigits()) return rc; | |
20 | ||
21 | return rc; | |
22 | ||
23 | } | |
24 | ||
25 | //_____________________________________________________________________________ | |
26 | Int_t AliTRDanalyzeHits() | |
27 | { | |
28 | // | |
29 | // Analyzes the hits and fills QA-histograms | |
30 | // | |
31 | ||
32 | Int_t rc = 0; | |
33 | ||
34 | //AliRunLoader *rl = gAlice->GetRunLoader(); | |
35 | AliRunLoader *rl = AliRunLoader::Open("TRD_test.root" | |
36 | ,AliConfig::GetDefaultEventFolderName()); | |
37 | if (!rl) { | |
38 | cout << "<AliTRDanalyzeHits> No RunLoader found" << endl; | |
39 | rc = 1; | |
40 | return rc; | |
41 | } | |
42 | ||
43 | AliLoader* loader = rl->GetLoader("TRDLoader"); | |
44 | if (!loader) { | |
45 | cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl; | |
46 | rc = 2; | |
47 | return rc; | |
48 | } | |
49 | ||
50 | rl->LoadgAlice(); | |
51 | rl->LoadHeader(); | |
52 | rl->LoadKinematics(); | |
53 | rl->LoadHits(); | |
54 | ||
55 | // Get the pointer to the TRD detector | |
56 | gAlice = rl->GetAliRun(); | |
57 | AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD"); | |
58 | if (!trd) { | |
59 | cout << "<AliTRDanalyzeHits> No TRD detector found" << endl; | |
60 | rc = 2; | |
61 | return rc; | |
62 | } | |
63 | ||
64 | // Get the pointer to the geometry object | |
65 | AliTRDgeometry *geo; | |
66 | if (trd) { | |
67 | geo = (AliTRDgeometry *) trd->GetGeometry(); | |
68 | } | |
69 | else { | |
70 | cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl; | |
71 | rc = 3; | |
72 | return rc; | |
73 | } | |
74 | ||
75 | AliTRDCommonParam *par = AliTRDCommonParam::Instance(); | |
76 | ||
77 | // Define the histograms | |
78 | TH1F *hQdedx = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0); | |
79 | TH1F *hQtr = new TH1F("hQtr" ,"Charge TR-hits" ,100,0.0,1000.0); | |
80 | ||
81 | Float_t rmin = geo->Rmin(); | |
82 | Float_t rmax = geo->Rmax(); | |
83 | Float_t length = geo->GetChamberLength(0,2); | |
84 | Float_t width = geo->GetChamberWidth(0); | |
85 | Int_t ncol = par->GetColMax(0); | |
86 | Int_t nrow = par->GetRowMax(0,2,13); | |
87 | Int_t ntime = ((Int_t) (rmax - rmin) / 22.0); | |
88 | ||
89 | TH2F *hZY = new TH2F("hZY" ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2. | |
90 | ,ntime, rmin, rmax); | |
91 | TH2F *hXZ = new TH2F("hXZ" ,"Z vs X (plane 0)" , ncol, -width/2., width/2. | |
92 | , nrow,-length/2.,length/2.); | |
93 | ||
94 | // Get the pointer hit tree | |
95 | TTree *hitTree = loader->TreeH(); | |
96 | if (!hitTree) { | |
97 | cout << "<AliTRDanalyzeHits> No hit tree found" << endl; | |
98 | rc = 4; | |
99 | return rc; | |
100 | } | |
101 | ||
102 | Int_t countHits = 0; | |
103 | Int_t nBytes = 0; | |
104 | ||
105 | // Get the number of entries in the hit tree | |
106 | // (Number of primary particles creating a hit somewhere) | |
107 | Int_t nTrack = (Int_t) hitTree->GetEntries(); | |
108 | cout << "<AliTRDanalyzeHits> Found " << nTrack | |
109 | << " primary particles with hits" << endl; | |
110 | ||
111 | // Loop through all entries in the tree | |
112 | for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) { | |
113 | ||
114 | gAlice->ResetHits(); | |
115 | nBytes += hitTree->GetEvent(iTrack); | |
116 | ||
117 | // Loop through the TRD hits | |
118 | Int_t iHit = 0; | |
119 | AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1); | |
120 | while (hit) { | |
121 | ||
122 | countHits++; | |
123 | iHit++; | |
124 | ||
125 | Float_t x = hit->X(); | |
126 | Float_t y = hit->Y(); | |
127 | Float_t z = hit->Z(); | |
128 | Float_t q = hit->GetCharge(); | |
129 | Int_t track = hit->Track(); | |
130 | Int_t det = hit->GetDetector(); | |
131 | Int_t plane = geo->GetPlane(det); | |
132 | ||
133 | if (q > 0) { | |
134 | hQdedx->Fill(q); | |
135 | hZY->Fill(z,y); | |
136 | if (plane == 0) { | |
137 | hXZ->Fill(x,z); | |
138 | } | |
139 | } | |
140 | else if (q < 0) { | |
141 | hQtr->Fill(TMath::Abs(q)); | |
142 | } | |
143 | ||
144 | hit = (AliTRDhit *) trd->NextHit(); | |
145 | ||
146 | } | |
147 | ||
148 | } | |
149 | ||
150 | cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl; | |
151 | ||
152 | TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600); | |
153 | cHits->Divide(2,2); | |
154 | cHits->cd(1); | |
155 | hXZ->Draw("COL"); | |
156 | cHits->cd(2); | |
157 | hZY->Draw("COL"); | |
158 | cHits->cd(3); | |
159 | gPad->SetLogy(); | |
160 | hQdedx->Draw(); | |
161 | cHits->cd(4); | |
162 | gPad->SetLogy(); | |
163 | hQtr->Draw(); | |
164 | ||
165 | return rc; | |
166 | ||
167 | } | |
168 | ||
169 | //_____________________________________________________________________________ | |
170 | Int_t AliTRDanalyzeDigits() | |
171 | { | |
172 | // | |
173 | // Analyzes the digits | |
174 | // | |
175 | ||
176 | Int_t rc = 0; | |
177 | ||
178 | const Int_t kNpla = AliTRDgeometry::Nplan(); | |
179 | ||
180 | if (!gAlice) { | |
181 | cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl; | |
182 | rc = 1; | |
183 | return rc; | |
184 | } | |
185 | ||
186 | AliRunLoader *rl = gAlice->GetRunLoader(); | |
187 | if (!rl) { | |
188 | cout << "<AliTRDanalyzeHits> No RunLoader found" << endl; | |
189 | rc = 2; | |
190 | return rc; | |
191 | } | |
192 | ||
193 | // Import the Trees for the event nEvent in the file | |
194 | rl->LoadDigits(); | |
195 | ||
196 | AliLoader* loader = rl->GetLoader("TRDLoader"); | |
197 | if (!loader) { | |
198 | cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl; | |
199 | rc = 3; | |
200 | return rc; | |
201 | } | |
202 | ||
203 | // Get the pointer to the TRD detector | |
204 | AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD"); | |
205 | if (!trd) { | |
206 | cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl; | |
207 | rc = 4; | |
208 | return rc; | |
209 | } | |
210 | ||
211 | // Get the parameter object | |
212 | AliTRDSimParam *parSim = AliTRDSimParam::Instance(); | |
213 | AliTRDCommonParam *parCom = AliTRDCommonParam::Instance(); | |
214 | AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); | |
215 | ||
216 | // Define the histograms | |
217 | Int_t adcRange = ((Int_t) parSim->GetADCoutRange()); | |
218 | TH1F *hAmpAll = new TH1F("hAmpAll" ,"Amplitude of the digits (all)" | |
219 | ,adcRange+1,-0.5,((Float_t) adcRange)+0.5); | |
220 | TH1F *hAmpEl = new TH1F("hAmpEl" ,"Amplitude of the digits (electrons)" | |
221 | ,adcRange+1,-0.5,((Float_t) adcRange)+0.5); | |
222 | TH1F *hAmpPi = new TH1F("hAmpPi" ,"Amplitude of the digits (pions)" | |
223 | ,adcRange+1,-0.5,((Float_t) adcRange)+0.5); | |
224 | TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)" | |
225 | ,5,-0.5,4.5); | |
226 | ||
227 | // Get the pointer to the geometry object | |
228 | AliTRDgeometry *geo; | |
229 | if (trd) { | |
230 | geo = trd->GetGeometry(); | |
231 | } | |
232 | else { | |
233 | cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl; | |
234 | rc = 5; | |
235 | return rc; | |
236 | } | |
237 | ||
238 | // Create the digits manager | |
239 | AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(); | |
240 | digitsManager->SetDebug(1); | |
241 | ||
242 | // Read the digits from the file | |
243 | if (!(digitsManager->ReadDigits(loader->TreeD()))) { | |
244 | cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl; | |
245 | rc = 6; | |
246 | return rc; | |
247 | } | |
248 | ||
249 | // Get the particle stack | |
250 | AliStack *kineStack = rl->Stack(); | |
251 | if (!kineStack) { | |
252 | cout << "<AliTRDanalyzeDigits> Cannot find the KINE stack" << endl; | |
253 | rc = 7; | |
254 | return rc; | |
255 | } | |
256 | ||
257 | Int_t countDigits = 0; | |
258 | Int_t iSec = 0; | |
259 | Int_t iCha = 2; | |
260 | Int_t timeMax = cal->GetNumberOfTimeBins(); | |
261 | ||
262 | TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)" | |
263 | ,timeMax,-0.5,((Double_t) timeMax)-0.5); | |
264 | TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)" | |
265 | ,timeMax,-0.5,((Double_t) timeMax)-0.5); | |
266 | ||
267 | // Loop over all planes | |
268 | for (Int_t iPla = 0; iPla < kNpla; iPla++) { | |
269 | ||
270 | Int_t iDet = geo->GetDetector(iPla,iCha,iSec); | |
271 | Int_t rowMax = parCom->GetRowMax(iPla,iCha,iSec); | |
272 | Int_t colMax = parCom->GetColMax(iPla); | |
273 | ||
274 | // Loop through the detector pixel | |
275 | for (Int_t time = 0; time < timeMax; time++) { | |
276 | for (Int_t col = 0; col < colMax; col++) { | |
277 | for (Int_t row = 0; row < rowMax; row++) { | |
278 | ||
279 | AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet); | |
280 | Int_t amp = digit->GetAmp(); | |
281 | Int_t track0 = digitsManager->GetTrack(0,row,col,time,iDet); | |
282 | Int_t track1 = digitsManager->GetTrack(1,row,col,time,iDet); | |
283 | TParticle *particle = 0; | |
284 | if (track0 > -1) { | |
285 | particle = (TParticle *) kineStack->Particle(track0); | |
286 | } | |
287 | ||
288 | if (amp > 0) { | |
289 | countDigits++; | |
290 | } | |
291 | ||
292 | // Total spectrum | |
293 | hAmpAll->Fill(amp); | |
294 | ||
295 | // Noise spectrum | |
296 | if (track0 < 0) { | |
297 | hAmpNoise->Fill(amp); | |
298 | } | |
299 | ||
300 | // Electron digit | |
301 | if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) { | |
302 | hAmpEl->Fill(amp); | |
303 | hAmpTimeEl->Fill(time,amp); | |
304 | } | |
305 | ||
306 | // Pion digit | |
307 | if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) { | |
308 | hAmpPi->Fill(amp); | |
309 | hAmpTimePi->Fill(time,amp); | |
310 | } | |
311 | ||
312 | delete digit; | |
313 | ||
314 | } | |
315 | } | |
316 | } | |
317 | ||
318 | } | |
319 | ||
320 | cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl; | |
321 | ||
322 | TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",100,100,600,800); | |
323 | cDigits->Divide(2,3); | |
324 | cDigits->cd(1); | |
325 | gPad->SetLogy(); | |
326 | hAmpAll->SetXTitle("Amplitude (ADC-channels)"); | |
327 | hAmpAll->SetYTitle("Entries"); | |
328 | hAmpAll->Draw(); | |
329 | cDigits->cd(2); | |
330 | gPad->SetLogy(); | |
331 | hAmpNoise->SetXTitle("Amplitude (ADC-channels)"); | |
332 | hAmpNoise->SetYTitle("Entries"); | |
333 | hAmpNoise->Draw(); | |
334 | cDigits->cd(3); | |
335 | gPad->SetLogy(); | |
336 | hAmpEl->SetXTitle("Amplitude (ADC-channels)"); | |
337 | hAmpEl->SetYTitle("Entries"); | |
338 | hAmpEl->Draw(); | |
339 | cDigits->cd(4); | |
340 | gPad->SetLogy(); | |
341 | hAmpPi->SetXTitle("Amplitude (ADC-channels)"); | |
342 | hAmpPi->SetYTitle("Entries"); | |
343 | hAmpPi->Draw(); | |
344 | cDigits->cd(5); | |
345 | hAmpTimeEl->SetXTitle("Timebin number"); | |
346 | hAmpTimeEl->SetYTitle("Mean amplitude"); | |
347 | hAmpTimeEl->Draw("HIST"); | |
348 | cDigits->cd(6); | |
349 | hAmpTimePi->SetXTitle("Timebin number"); | |
350 | hAmpTimePi->SetYTitle("Mean amplitude"); | |
351 | hAmpTimePi->Draw("HIST"); | |
352 | ||
353 | return rc; | |
354 | ||
355 | } | |
356 |