3 void MUONdoubles (Int_t evNumber1=0,Int_t evNumber2=0)
6 hprof = new TProfile("hprof","Profile dmin vs r",24,0,120,0,50);
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++) {
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));
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;
29 for (Int_t i=0; i<10; i++) {
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]));
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.
48 /////////////////////////////////////////////////////////////////////////
50 // Dynamically link some shared libs
52 if (gClassTable->GetID("AliRun") < 0) {
53 gROOT->LoadMacro("loadlibs.C");
57 // Connect the Root Galice file containing Geometry, Kine and Hits
59 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
60 if (!file) file = new TFile("galice.root","UPDATE");
62 // Get AliRun object from file or create it if not on file
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");
70 // Set reconstruction models
72 // Get pointers to Alice detectors and Digits containers
73 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
75 // Book some histograms
76 TH1F *Hcentre[10], *Hall[10], *Hcfrac[10], *Htail[10], *Htfrac[10];
77 TH1F *Hdcut[10], *Hdall[10], *Hdfrac[10];
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)");
84 Hall[i] = new TH1F("Hall","Hit Distribution",26,0,260);
85 Hall[i]->SetFillColor(0);
86 Hall[i]->SetXTitle("R (cm)");
88 Hcfrac[i] = new TH1F("Hcfrag","Hit Distribution",26,0,260);
89 Hcfrac[i]->SetFillColor(0);
90 Hcfrac[i]->SetXTitle("R (cm)");
92 Htail[i] = new TH1F("Htail","Hit Distribution",26,0,260);
93 Htail[i]->SetFillColor(0);
94 Htail[i]->SetXTitle("R (cm)");
96 Htfrac[i] = new TH1F("Htfrag","Hit Distribution",26,0,260);
97 Htfrac[i]->SetFillColor(0);
98 Htfrac[i]->SetXTitle("R (cm)");
100 Hdall[i] = new TH1F("Hdall","Hit Distribution",26,0,260);
101 Hdall[i]->SetFillColor(0);
102 Hdall[i]->SetXTitle("R (cm)");
104 Hdcut[i] = new TH1F("Hdcut","Hit Distribution",26,0,260);
105 Hdcut[i]->SetFillColor(0);
106 Hdcut[i]->SetXTitle("R (cm)");
108 Hdfrac[i] = new TH1F("Hdfrac","Hit Distribution",26,0,260);
109 Hdfrac[i]->SetFillColor(0);
110 Hdfrac[i]->SetXTitle("R (cm)");
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;
126 TTree *TH = gAlice->TreeH();
127 Int_t ntracks = TH->GetEntries();
128 cout<<"ntracks "<<ntracks<<endl;
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);
143 // Loop on chambers and on cathode planes
145 TTree *TH = gAlice->TreeH();
146 Int_t ntracks = TH->GetEntries();
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.;
160 for (Int_t track=0; track<ntracks; track++) {
162 Int_t nbytes += TH->GetEvent(track);
163 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
165 mHit=(AliMUONHit*)MUON->NextHit())
167 Int_t nch = mHit->fChamber; // chamber number
168 if (nch!=1) continue;
170 x[nhit] = mHit->fX; // x-pos of hit
171 y[nhit] = mHit->fY; // y-pos
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;
183 if (dx<0.5 && dy < 0.75) Hdcut[0]->Fill(r,1.);
186 if (r>20) hprof->Fill(r,dmin,1);
189 gAlice->ResetDigits();
190 gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
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;
197 // Get ready the current chamber stuff
199 AliMUONResponse* response = iChamber.GetResponseModel();
200 AliMUONSegmentation* seg = iChamber.GetSegmentationModel(icat+1);
201 HitMap = new AliMUONHitMapA1(seg, MUONdigits);
203 Int_t nxmax=seg->Npx();
204 Int_t nymax=seg->Npy();
206 // generate random positions on the chamber
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);
212 // test for centre overlap
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.);
222 // test for neighbour overlap
226 seg->Neighbours(ix,iy,&nn,X,Y);
229 for (Int_t j=0; j<nn; j++) {
230 if (HitMap->TestHit(X[j],Y[j]) != 0) hit=true;
232 if (hit && HitMap->TestHit(ix,iy) == 0) Htail[ich]->Fill(r,1.);
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]);
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);
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);