]>
Commit | Line | Data |
---|---|---|
ddae0931 | 1 | void RICHtest (Int_t evNumber1=0,Int_t evNumber2=0) |
2 | { | |
3 | ///////////////////////////////////////////////////////////////////////// | |
4 | // This macro is a small example of a ROOT macro | |
5 | // illustrating how to read the output of GALICE | |
6 | // and do some analysis. | |
7 | // | |
8 | ///////////////////////////////////////////////////////////////////////// | |
9 | ||
10 | // Dynamically link some shared libs | |
11 | ||
12 | if (gClassTable->GetID("AliRun") < 0) { | |
13 | gROOT->LoadMacro("loadlibs.C"); | |
14 | loadlibs(); | |
15 | } | |
16 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
17 | ||
18 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
19 | ||
20 | ||
21 | if (!file) { | |
22 | printf("\n Creating galice.root \n"); | |
23 | file = new TFile("galice.root"); | |
24 | } else { | |
25 | printf("\n galice.root found in file list"); | |
26 | } | |
27 | file->ls(); | |
28 | ||
29 | // Get AliRun object from file or create it if not on file | |
30 | if (!gAlice) { | |
31 | gAlice = (AliRun*)(file->Get("gAlice")); | |
32 | if (gAlice) printf("AliRun object found on file\n"); | |
33 | if (!gAlice) { | |
34 | printf("\n create new gAlice object"); | |
35 | gAlice = new AliRun("gAlice","Alice test program"); | |
36 | } | |
37 | } | |
38 | ||
39 | // Create some histograms | |
40 | ||
41 | ||
42 | TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution" | |
43 | ,100,0.,500.); | |
44 | TH1F *pcharge1 = new TH1F("pcharge1","Pad Charge distribution" | |
45 | ,100,0.,200.); | |
46 | TH1F *xresid1 = new TH1F("xresid1","x-Residuals" | |
47 | ,100,-0.5,0.5); | |
48 | TH1F *yresid1 = new TH1F("yresid1","y-Residuals" | |
49 | ,100,-0.2,0.2); | |
50 | TH1F *npads1 = new TH1F("npads1" ,"Pads per Hit" | |
51 | ,20,-0.5,19.5); | |
52 | TH2F *xresys1 = new TH2F("xresys1","x-Residuals systematics" | |
53 | ,50,-2.0,2.0,100,-1.0,1.0); | |
54 | TH2F *yresys1 = new TH2F("yresys1","y-Residuals systematics" | |
55 | ,50,-0.4,0.4,100,-0.1,0.1); | |
56 | ||
57 | TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution" | |
58 | ,100,0.,500.); | |
59 | TH1F *pcharge2 = new TH1F("pcharge2","Pad Charge distribution" | |
60 | ,100,0.,200.); | |
61 | TH1F *xresid2 = new TH1F("xresid2","x-Residuals" | |
62 | ,100,-0.5,0.5); | |
63 | TH1F *yresid2 = new TH1F("yresid2","y-Residuals" | |
64 | ,100,-0.2,0.2); | |
65 | TH1F *npads2 = new TH1F("npads2" ,"Pads per Hit" | |
66 | ,20,-0.5,19.5); | |
67 | TH2F *xresys2 = new TH2F("xresys2","x-Residuals systematics" | |
68 | ,50,-2.0,2.0,100,-1.0,1.0); | |
69 | TH2F *yresys2 = new TH2F("yresys2","y-Residuals systematics" | |
70 | ,50,-0.4,0.4,100,-0.1,0.1); | |
71 | ||
72 | ||
73 | AliRICHchamber* iChamber; | |
74 | // | |
75 | // Loop over events | |
76 | // | |
77 | Int_t Nh=0; | |
78 | Int_t Nh1=0; | |
79 | for (Int_t nev=0; nev<= evNumber2; nev++) { | |
80 | //cout << "nev " << nev <<endl; | |
81 | printf ("Event Number : %d\n",nev); | |
82 | Int_t nparticles = gAlice->GetEvent(nev); | |
83 | //cout << "nparticles " << nparticles <<endl; | |
84 | printf ("Number of particles: %d\n",nparticles); | |
85 | if (nev < evNumber1) continue; | |
86 | if (nparticles <= 0) return; | |
87 | ||
88 | AliRICH *RICH = (AliRICH*) gAlice->GetDetector("RICH"); | |
89 | printf("\n track %d %d \n ", nev, RICH); | |
90 | ||
91 | TTree *TH = gAlice->TreeH(); | |
92 | Int_t ntracks = TH->GetEntries(); | |
93 | Int_t Nc=0; | |
94 | Int_t npad[2]; | |
95 | Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2]; | |
96 | ||
97 | printf ("Just got it"); | |
98 | // | |
99 | // Loop over events | |
100 | // | |
101 | for (Int_t track=0; track<ntracks;track++) { | |
102 | ||
103 | gAlice->ResetHits(); | |
104 | Int_t nbytes += TH->GetEvent(track); | |
105 | if (RICH) { | |
106 | for(AliRICHhit* mHit=(AliRICHhit*)RICH->FirstHit(-1); | |
107 | mHit; | |
108 | mHit=(AliRICHhit*)RICH->NextHit()) | |
109 | { | |
110 | Int_t nch = mHit->fChamber; // chamber number | |
111 | Float_t x = mHit->fX; // x-pos of hit | |
112 | Float_t y = mHit->fY; // y-pos | |
113 | // | |
114 | // | |
115 | iChamber = & RICH->Chamber(nch-1); | |
116 | response=iChamber->GetResponseModel(mip); | |
117 | // | |
118 | // | |
119 | if (nch <= 1) { | |
120 | for (Int_t i=0; i<2; i++) { | |
121 | xmean[i]=0; | |
122 | ymean[i]=0; | |
123 | Q[i]=0; | |
124 | npad[i]=0; | |
125 | } | |
126 | for (AliRICHcluster *mPad=(AliRICHcluster*)RICH | |
127 | ->FirstPad(mHit); | |
128 | mPad; | |
129 | mPad=(AliRICHcluster*)RICH->NextPad() | |
130 | ) | |
131 | { | |
132 | Int_t nseg = mPad->fCathode; // segmentation | |
133 | Int_t nhit = mPad->fHitNumber; // hit number | |
134 | Int_t qtot = mPad->fQ; // charge | |
135 | Int_t ipx = mPad->fPadX; // pad number on X | |
136 | Int_t ipy = mPad->fPadY; // pad number on Y | |
137 | Int_t iqpad = mPad->fQpad; // charge per pad | |
138 | Int_t secpad = mPad->fRSec; // r-pos of pad | |
139 | // | |
140 | // | |
141 | segmentation=iChamber->GetSegmentationModel(nseg); | |
142 | Int_t ind=nseg-1; | |
143 | // printf("Pad of hit, padx, pady, iqpad, ncha %d %d %d %d %d %d\n", | |
144 | // nhit, ipx, ipy, iqpad, nseg, ind); | |
145 | ||
146 | ||
147 | if (nseg == 1) { | |
148 | pcharge1->Fill(iqpad,(float) 1); | |
149 | ccharge1->Fill(qtot ,(float) 1); | |
150 | } else { | |
151 | pcharge2->Fill(iqpad,(float) 1); | |
152 | ccharge2->Fill(qtot ,(float) 1); | |
153 | } | |
154 | // Calculate centre of gravity | |
155 | // | |
156 | if (iqpad == 0 && ind==0) { | |
157 | printf("\n attention iqpad 0"); | |
158 | } | |
159 | ||
160 | if (iqpad > 0 && secpad==1) { | |
161 | Float_t xc,yc; | |
162 | npad[ind]++; | |
163 | segmentation->GetPadCxy(ipx,ipy,xc,yc); | |
164 | // printf("\n pad %d %d ", ipx, ipy); | |
165 | // printf("\n pos %f %f ", xc, yc); | |
166 | xmean[ind]+=Float_t(iqpad*xc); | |
167 | ymean[ind]+=Float_t(iqpad*yc); | |
168 | Q[ind]+=iqpad; | |
169 | } | |
170 | ||
171 | } //pad hit loop | |
172 | for (Int_t i=0; i<2; i++) { | |
173 | segmentation = iChamber->GetSegmentationModel(i+1); | |
174 | ||
175 | // if (npad[i] ==0) { | |
176 | // printf("\n #pads=0 %f %f",x,y); | |
177 | // } | |
178 | ||
179 | if (Q[i] >0) { | |
180 | xmean[i] = xmean[i]/Q[i]; | |
181 | xres[i] = xmean[i]-x; | |
182 | ymean[i] = ymean[i]/Q[i]; | |
183 | yres[i] = ymean[i]-y; | |
184 | // printf("\n xmean, ymean %f %f",xmean[i],ymean[i]); | |
185 | // printf("\n xhit, yhit %f %f",x,y); | |
186 | ||
187 | // Systematic Error | |
188 | // | |
189 | Int_t icx, icy; | |
190 | segmentation-> | |
191 | GetPadIxy(x,y,icx,icy); | |
192 | segmentation-> | |
193 | GetPadCxy(icx,icy,xonpad[i],yonpad[i]); | |
194 | xonpad[i]-=x; | |
195 | yonpad[i]-=y; | |
196 | } | |
197 | } // plane loop | |
198 | xresid1->Fill(xres[0],(float) 1); | |
199 | yresid1->Fill(yres[0],(float) 1); | |
200 | npads1->Fill(npad[0],(float) 1); | |
201 | if (npad[0] >=2 && Q[0] > 6 ) { | |
202 | xresys1->Fill(xonpad[0],xres[0],(float) 1); | |
203 | yresys1->Fill(yonpad[0],yres[0],(float) 1); | |
204 | } | |
205 | ||
206 | xresid2->Fill(xres[1],(float) 1); | |
207 | yresid2->Fill(yres[1],(float) 1); | |
208 | npads2->Fill(npad[1],(float) 1); | |
209 | if (npad[1] >=2 && Q[1] > 6) { | |
210 | xresys2->Fill(xonpad[1],xres[1],(float) 1); | |
211 | yresys2->Fill(yonpad[1],yres[1],(float) 1); | |
212 | } | |
213 | } // chamber 1 | |
214 | } // hit loop | |
215 | } // if RICH | |
216 | } // track loop | |
217 | } // event loop | |
218 | ||
219 | //Create a canvas, set the view range, show histograms | |
220 | TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700); | |
221 | pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99); | |
222 | pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99); | |
223 | pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49); | |
224 | pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49); | |
225 | pad11->SetFillColor(11); | |
226 | pad12->SetFillColor(11); | |
227 | pad13->SetFillColor(11); | |
228 | pad14->SetFillColor(11); | |
229 | pad11->Draw(); | |
230 | pad12->Draw(); | |
231 | pad13->Draw(); | |
232 | pad14->Draw(); | |
233 | ||
234 | pad11->cd(); | |
235 | ccharge1->SetFillColor(42); | |
236 | ccharge1->SetXTitle("ADC units"); | |
237 | ccharge1->Draw(); | |
238 | ||
239 | pad12->cd(); | |
240 | pcharge1->SetFillColor(42); | |
241 | pcharge1->SetXTitle("ADC units"); | |
242 | pcharge1->Draw(); | |
243 | ||
244 | pad13->cd(); | |
245 | xresid1->SetFillColor(42); | |
246 | xresid1->Draw(); | |
247 | ||
248 | pad14->cd(); | |
249 | yresid1->SetFillColor(42); | |
250 | yresid1->Draw(); | |
251 | ||
252 | TCanvas *c2 = new TCanvas("c2","Cluster Size",400,10,600,700); | |
253 | pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99); | |
254 | pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99); | |
255 | pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49); | |
256 | pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49); | |
257 | pad21->SetFillColor(11); | |
258 | pad22->SetFillColor(11); | |
259 | pad23->SetFillColor(11); | |
260 | pad24->SetFillColor(11); | |
261 | pad21->Draw(); | |
262 | pad22->Draw(); | |
263 | pad23->Draw(); | |
264 | pad24->Draw(); | |
265 | ||
266 | pad21->cd(); | |
267 | npads1->SetFillColor(42); | |
268 | npads1->SetXTitle("Cluster Size"); | |
269 | npads1->Draw(); | |
270 | ||
271 | pad23->cd(); | |
272 | xresys1->SetXTitle("x on pad"); | |
273 | xresys1->SetYTitle("x-xcog"); | |
274 | xresys1->Draw(); | |
275 | ||
276 | pad24->cd(); | |
277 | yresys1->SetXTitle("y on pad"); | |
278 | yresys1->SetYTitle("y-ycog"); | |
279 | yresys1->Draw(); | |
280 | ||
281 | TCanvas *c3 = new TCanvas("c3","Charge and Residuals",400,10,600,700); | |
282 | pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99); | |
283 | pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99); | |
284 | pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49); | |
285 | pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49); | |
286 | pad31->SetFillColor(11); | |
287 | pad32->SetFillColor(11); | |
288 | pad33->SetFillColor(11); | |
289 | pad34->SetFillColor(11); | |
290 | pad31->Draw(); | |
291 | pad32->Draw(); | |
292 | pad33->Draw(); | |
293 | pad34->Draw(); | |
294 | ||
295 | pad31->cd(); | |
296 | ccharge2->SetFillColor(42); | |
297 | ccharge2->SetXTitle("ADC units"); | |
298 | ccharge2->Draw(); | |
299 | ||
300 | pad32->cd(); | |
301 | pcharge2->SetFillColor(42); | |
302 | pcharge2->SetXTitle("ADC units"); | |
303 | pcharge2->Draw(); | |
304 | ||
305 | pad33->cd(); | |
306 | xresid2->SetFillColor(42); | |
307 | xresid2->Draw(); | |
308 | ||
309 | pad34->cd(); | |
310 | yresid2->SetFillColor(42); | |
311 | yresid2->Draw(); | |
312 | ||
313 | TCanvas *c4 = new TCanvas("c4","Cluster Size",400,10,600,700); | |
314 | pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99); | |
315 | pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99); | |
316 | pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49); | |
317 | pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49); | |
318 | pad41->SetFillColor(11); | |
319 | pad42->SetFillColor(11); | |
320 | pad43->SetFillColor(11); | |
321 | pad44->SetFillColor(11); | |
322 | pad41->Draw(); | |
323 | pad42->Draw(); | |
324 | pad43->Draw(); | |
325 | pad44->Draw(); | |
326 | ||
327 | pad41->cd(); | |
328 | npads2->SetFillColor(42); | |
329 | npads2->SetXTitle("Cluster Size"); | |
330 | npads2->Draw(); | |
331 | ||
332 | pad43->cd(); | |
333 | xresys2->SetXTitle("x on pad"); | |
334 | xresys2->SetYTitle("x-xcog"); | |
335 | xresys2->Draw(); | |
336 | ||
337 | pad44->cd(); | |
338 | yresys2->SetXTitle("y on pad"); | |
339 | yresys2->SetYTitle("y-ycog"); | |
340 | yresys2->Draw(); | |
341 | ||
342 | } |