]>
Commit | Line | Data |
---|---|---|
a9e2aefa | 1 | #include "iostream.h" |
2 | ||
3 | void MUONdoubles (Int_t evNumber1=0,Int_t evNumber2=0) | |
4 | { | |
5 | { | |
6 | hprof = new TProfile("hprof","Profile dmin vs r",24,0,120,0,50); | |
7 | } | |
8 | hdist1 = new TH1F("hdist1","distance",100,0,10); | |
9 | hdist2 = new TH1F("hdist2","distance",100,0,10); | |
10 | Float_t a1=3.7/TMath::Sqrt(4); | |
11 | Float_t a2=26/TMath::Sqrt(4); | |
12 | for (Int_t i=1; i<10000; i++) { | |
13 | Float_t x1,x2,y1,y2; | |
14 | x1 = (2*gRandom->Rndm(i)-1.)*a1; | |
15 | x2 = (2*gRandom->Rndm(i)-1.)*a1; | |
16 | y1 = (2*gRandom->Rndm(i)-1.)*a1; | |
17 | y2 = (2*gRandom->Rndm(i)-1.)*a1; | |
18 | Float_t d=TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); | |
19 | hdist1->Fill(d,1); | |
20 | } | |
21 | Float_t xh[100], yh[100]; | |
22 | for (Int_t k=0; k<1000; k++) { | |
23 | for (Int_t i=0; i<100; i++) { | |
24 | xh[i] = (2*gRandom->Rndm(i)-1.)*a2; | |
25 | yh[i] = (2*gRandom->Rndm(i)-1.)*a2; | |
26 | } | |
27 | Float_t x1,y1; | |
28 | ||
29 | for (Int_t i=0; i<10; i++) { | |
30 | Float_t dmin=1000; | |
31 | x1 = (2*gRandom->Rndm(i)-1.)*a2; | |
32 | y1 = (2*gRandom->Rndm(i)-1.)*a2; | |
33 | for (Int_t j=0; j<100; j++) { | |
34 | Float_t d=TMath::Sqrt((x1-xh[j])*(x1-xh[j])+(y1-yh[j])*(y1-yh[j])); | |
35 | if (d<dmin) dmin=d; | |
36 | } | |
37 | hdist2->Fill(dmin,1); | |
38 | } | |
39 | } | |
40 | ||
41 | ||
42 | ||
43 | ///////////////////////////////////////////////////////////////////////// | |
44 | // This macro is a small example of a ROOT macro | |
45 | // illustrating how to read the output of GALICE | |
46 | // and do some analysis. | |
47 | // | |
48 | ///////////////////////////////////////////////////////////////////////// | |
49 | ||
50 | // Dynamically link some shared libs | |
51 | ||
52 | if (gClassTable->GetID("AliRun") < 0) { | |
53 | gROOT->LoadMacro("loadlibs.C"); | |
54 | loadlibs(); | |
55 | } | |
56 | ||
57 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
58 | ||
59 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
60 | if (!file) file = new TFile("galice.root","UPDATE"); | |
61 | ||
62 | // Get AliRun object from file or create it if not on file | |
63 | ||
64 | if (!gAlice) { | |
65 | gAlice = (AliRun*)file->Get("gAlice"); | |
66 | if (gAlice) printf("AliRun object found on file\n"); | |
67 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
68 | } | |
69 | // | |
70 | // Set reconstruction models | |
71 | // | |
72 | // Get pointers to Alice detectors and Digits containers | |
73 | AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); | |
74 | // | |
75 | // Book some histograms | |
76 | TH1F *Hcentre[10], *Hall[10], *Hcfrac[10], *Htail[10], *Htfrac[10]; | |
77 | TH1F *Hdcut[10], *Hdall[10], *Hdfrac[10]; | |
78 | ||
79 | for (Int_t i=0; i<10; i++) { | |
80 | Hcentre[i] = new TH1F("Hcentre","Hit Distribution",26,0,260); | |
81 | Hcentre[i]->SetFillColor(0); | |
82 | Hcentre[i]->SetXTitle("R (cm)"); | |
83 | ||
84 | Hall[i] = new TH1F("Hall","Hit Distribution",26,0,260); | |
85 | Hall[i]->SetFillColor(0); | |
86 | Hall[i]->SetXTitle("R (cm)"); | |
87 | ||
88 | Hcfrac[i] = new TH1F("Hcfrag","Hit Distribution",26,0,260); | |
89 | Hcfrac[i]->SetFillColor(0); | |
90 | Hcfrac[i]->SetXTitle("R (cm)"); | |
91 | ||
92 | Htail[i] = new TH1F("Htail","Hit Distribution",26,0,260); | |
93 | Htail[i]->SetFillColor(0); | |
94 | Htail[i]->SetXTitle("R (cm)"); | |
95 | ||
96 | Htfrac[i] = new TH1F("Htfrag","Hit Distribution",26,0,260); | |
97 | Htfrac[i]->SetFillColor(0); | |
98 | Htfrac[i]->SetXTitle("R (cm)"); | |
99 | ||
100 | Hdall[i] = new TH1F("Hdall","Hit Distribution",26,0,260); | |
101 | Hdall[i]->SetFillColor(0); | |
102 | Hdall[i]->SetXTitle("R (cm)"); | |
103 | ||
104 | Hdcut[i] = new TH1F("Hdcut","Hit Distribution",26,0,260); | |
105 | Hdcut[i]->SetFillColor(0); | |
106 | Hdcut[i]->SetXTitle("R (cm)"); | |
107 | ||
108 | Hdfrac[i] = new TH1F("Hdfrac","Hit Distribution",26,0,260); | |
109 | Hdfrac[i]->SetFillColor(0); | |
110 | Hdfrac[i]->SetXTitle("R (cm)"); | |
111 | ||
112 | } | |
113 | ||
114 | // | |
115 | // Loop over events | |
116 | // | |
117 | Int_t Nh=0; | |
118 | Int_t Nh1=0; | |
119 | for (int nev=0; nev<= evNumber2; nev++) { | |
120 | Int_t nparticles = gAlice->GetEvent(nev); | |
121 | cout << "nev " << nev <<endl; | |
122 | cout << "nparticles " << nparticles <<endl; | |
123 | if (nev < evNumber1) continue; | |
124 | if (nparticles <= 0) return; | |
125 | ||
126 | TTree *TH = gAlice->TreeH(); | |
127 | Int_t ntracks = TH->GetEntries(); | |
128 | cout<<"ntracks "<<ntracks<<endl; | |
129 | ||
130 | Int_t nbytes = 0; | |
131 | ||
132 | ||
133 | TClonesArray *Particles = gAlice->Particles(); | |
134 | TTree *TD = gAlice->TreeD(); | |
135 | Int_t nent=TD->GetEntries(); | |
136 | printf("Found %d entries in the tree (must be one per cathode per event!)\n",nent); | |
137 | Float_t x[2000]; | |
138 | Float_t y[2000]; | |
139 | Int_t nhit=0; | |
140 | ||
141 | if (MUON) { | |
142 | // | |
143 | // Loop on chambers and on cathode planes | |
144 | // | |
145 | TTree *TH = gAlice->TreeH(); | |
146 | Int_t ntracks = TH->GetEntries(); | |
147 | ||
148 | // | |
149 | // Loop over Hits | |
150 | // | |
151 | for (i=0; i<1000; i++) { | |
152 | Float_t r =100.*gRandom->Rndm(i); | |
153 | Float_t phi=2.*TMath::Pi()*gRandom->Rndm(i); | |
154 | Float_t xr=r*TMath::Sin(phi); | |
155 | Float_t yr=r*TMath::Cos(phi); | |
156 | Hdall[0]->Fill(r,1.); | |
157 | Float_t dmin=100000.; | |
158 | ||
159 | if (i==0) { | |
160 | for (Int_t track=0; track<ntracks; track++) { | |
161 | gAlice->ResetHits(); | |
162 | Int_t nbytes += TH->GetEvent(track); | |
163 | for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); | |
164 | mHit; | |
165 | mHit=(AliMUONHit*)MUON->NextHit()) | |
166 | { | |
167 | Int_t nch = mHit->fChamber; // chamber number | |
168 | if (nch!=1) continue; | |
169 | ||
170 | x[nhit] = mHit->fX; // x-pos of hit | |
171 | y[nhit] = mHit->fY; // y-pos | |
172 | nhit++; | |
173 | } // hit loop | |
174 | } // track loop | |
175 | } else { | |
176 | for (Int_t nh=0; nh<nhit; nh++) { | |
177 | Float_t d=TMath::Sqrt((x[nh]-xr)*(x[nh]-xr)+(y[nh]-yr)*(y[nh]-yr)); | |
178 | // printf ("\n r,d %f %f",r,d); | |
179 | Float_t dx=TMath::Abs(x[nh]-xr); | |
180 | Float_t dy=TMath::Abs(y[nh]-yr); | |
181 | if (d < dmin) dmin=d; | |
182 | ||
183 | if (dx<0.5 && dy < 0.75) Hdcut[0]->Fill(r,1.); | |
184 | } // hit loop | |
185 | } // i==0 | |
186 | if (r>20) hprof->Fill(r,dmin,1); | |
187 | } // random loop | |
188 | Int_t icat=0; | |
189 | gAlice->ResetDigits(); | |
190 | gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ... | |
191 | ||
192 | for (Int_t ich=0;ich<1;ich++) { | |
193 | AliMUONChamber iChamber= MUON->Chamber(ich); | |
194 | TClonesArray *MUONdigits = MUON->DigitsAddress(ich); | |
195 | if (MUONdigits == 0) continue; | |
196 | // | |
197 | // Get ready the current chamber stuff | |
198 | // | |
199 | AliMUONResponse* response = iChamber.GetResponseModel(); | |
200 | AliMUONSegmentation* seg = iChamber.GetSegmentationModel(icat+1); | |
201 | HitMap = new AliMUONHitMapA1(seg, MUONdigits); | |
202 | HitMap->FillHits(); | |
203 | Int_t nxmax=seg->Npx(); | |
204 | Int_t nymax=seg->Npy(); | |
205 | // | |
206 | // generate random positions on the chamber | |
207 | // | |
208 | for (Int_t i=0; i<4000; i++) { | |
209 | Int_t ix = 2*Int_t((gRandom->Rndm(i)-0.5)*nxmax); | |
210 | Int_t iy = 2*Int_t((gRandom->Rndm(i)-0.5)*nymax); | |
211 | ||
212 | // test for centre overlap | |
213 | // | |
214 | Float_t xp,yp; | |
215 | seg->GetPadCxy(ix,iy,xp,yp); | |
216 | Float_t r=TMath::Sqrt(xp*xp+yp*yp); | |
217 | if (TMath::Abs(xp)>3 && TMath::Abs(yp)>3) | |
218 | Hall[ich]->Fill(r,1.); | |
219 | if (HitMap->TestHit(ix,iy) != 0) { | |
220 | Hcentre[ich]->Fill(r,1.); | |
221 | } | |
222 | // test for neighbour overlap | |
223 | // | |
224 | Int_t nn; | |
225 | Int_t X[20], Y[20]; | |
226 | seg->Neighbours(ix,iy,&nn,X,Y); | |
227 | Bool_t hit=false; | |
228 | ||
229 | for (Int_t j=0; j<nn; j++) { | |
230 | if (HitMap->TestHit(X[j],Y[j]) != 0) hit=true; | |
231 | } | |
232 | if (hit && HitMap->TestHit(ix,iy) == 0) Htail[ich]->Fill(r,1.); | |
233 | } //random loop | |
234 | delete HitMap; | |
235 | } // chamber loop | |
236 | } // event loop | |
237 | } // if MUON | |
238 | ||
239 | for (Int_t ich=0; ich<10; ich++) { | |
240 | Hcfrac[ich]->Divide(Hcentre[ich],Hall[ich]); | |
241 | Htfrac[ich]->Divide(Htail[ich],Hall[ich]); | |
242 | Hdfrac[ich]->Divide(Hdcut[ich],Hdall[ich]); | |
243 | } | |
244 | ||
245 | TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700); | |
246 | pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99); | |
247 | pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99); | |
248 | pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49); | |
249 | pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49); | |
250 | pad11->SetFillColor(0); | |
251 | pad12->SetFillColor(0); | |
252 | pad13->SetFillColor(0); | |
253 | pad14->SetFillColor(0); | |
254 | pad11->Draw(); | |
255 | pad12->Draw(); | |
256 | pad13->Draw(); | |
257 | pad14->Draw(); | |
258 | ||
259 | pad11->cd(); | |
260 | Hcentre[0]->Draw(); | |
261 | ||
262 | pad12->cd(); | |
263 | Hall[0]->Draw(); | |
264 | ||
265 | pad13->cd(); | |
266 | Hcfrac[0]->Draw(); | |
267 | ||
268 | pad14->cd(); | |
269 | Htfrac[0]->Draw(); | |
270 | ||
271 | TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700); | |
272 | pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99); | |
273 | pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99); | |
274 | pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49); | |
275 | pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49); | |
276 | pad21->SetFillColor(0); | |
277 | pad22->SetFillColor(0); | |
278 | pad23->SetFillColor(0); | |
279 | pad24->SetFillColor(0); | |
280 | pad21->Draw(); | |
281 | pad22->Draw(); | |
282 | pad23->Draw(); | |
283 | pad24->Draw(); | |
284 | ||
285 | pad21->cd(); | |
286 | hdist1->Draw(); | |
287 | ||
288 | pad22->cd(); | |
289 | hdist2->Draw(); | |
290 | ||
291 | pad23->cd(); | |
292 | Hdfrac[0]->Draw(); | |
293 | ||
294 | pad24->cd(); | |
295 | hprof->Draw(); | |
296 | file->Close(); | |
297 | } | |
298 |