Code from MUON-dev joined
[u/mrichter/AliRoot.git] / MUON / MUONhitrectest.C
CommitLineData
a9e2aefa 1#include "iostream.h"
2
3void MUONhitrectest (Int_t evNumber1=0,Int_t evNumber2=0)
4{
5/////////////////////////////////////////////////////////////////////////
6// This macro is a small example of a ROOT macro
7// illustrating how to read the output of GALICE
8// and do some analysis.
9//
10/////////////////////////////////////////////////////////////////////////
11
12// Dynamically link some shared libs
13
14 if (gClassTable->GetID("AliRun") < 0) {
15 gROOT->LoadMacro("loadlibs.C");
16 loadlibs();
17 }
18
19// Connect the Root Galice file containing Geometry, Kine and Hits
20
21 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
22 if (!file) file = new TFile("galice.root");
23
24// Get AliRun object from file or create it if not on file
25
26 if (!gAlice) {
27 gAlice = (AliRun*)file->Get("gAlice");
28 if (gAlice) printf("AliRun object found on file\n");
29 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
30 }
31// Create some histograms
32
33 TH2F *h21 = new TH2F("h21","Hits",100,-100,100,100,-100,100);
34 TH2F *h22 = new TH2F("h22","CoG ",100,-100,100,100,-100,100);
35 TH1F *h1 = new TH1F("h1","Multiplicity",30,-0.5,29.5);
36 TH1F *hmult = new TH1F("hmult","Multiplicity",30,-0.5,29.5);
37 TH1F *hresx = new TH1F("hresx","Residuals",100,-1,1);
38 TH1F *hresy = new TH1F("hresy","Residuals",100,-10.,10);
39 TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
40 TH1F *hpos = new TH1F("hnpos","Possibilities",10,-0.5,9.5);
41 TH2F *hchi1 = new TH2F("hchi1","Chi2 vs Residuals",100,0,0.2,100,-500,500);
42 TH2F *hchi2 = new TH2F("hchi2","Chi2 vs Residuals",100,0,20,100,-500,500);
43//
44// Loop over events
45//
46 Int_t Nh=0;
47 Int_t Nh1=0;
48
49 for (int nev=0; nev<= evNumber2; nev++) {
50 Int_t nparticles = gAlice->GetEvent(nev);
51 cout << "nev " << nev <<endl;
52 cout << "nparticles " << nparticles <<endl;
53 if (nev < evNumber1) continue;
54 if (nparticles <= 0) return;
55
56 TTree *TH = gAlice->TreeH();
57 Int_t ntracks = TH->GetEntries();
58 cout<<"ntracks "<<ntracks<<endl;
59
60 Int_t nbytes = 0;
61
62
63
64// Get pointers to Alice detectors and Digits containers
65 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
66 TClonesArray *Particles = gAlice->Particles();
67 TTree *TR = gAlice->TreeR();
68 Int_t nent=TR->GetEntries();
69 printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
70 if (MUON) {
71 Float_t xc[2000], yc[2000], itrc[2000];
72 TClonesArray *MUONrawclust = MUON->RawClustAddress(0);
73 nbytes += TR->GetEvent(1);
74 Int_t nrawcl = MUONrawclust->GetEntries();
75 AliMUONRawCluster *mRaw;
76 for (Int_t iraw=0; iraw < nrawcl; iraw++) {
77 mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
78 xc[iraw]=mRaw->fX[1];
79 yc[iraw]=mRaw->fY[0];
80 itrc[iraw]=mRaw->fTracks[1];
81 } // cluster
82
83 gAlice->ResetHits();
84 for (Int_t itrack=0; itrack<ntracks ;itrack++)
85 {
86 Int_t nbytes += TH->GetEvent(itrack);
87 Int_t nhit=MUON->Hits()->GetEntriesFast();
88 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(itrack);
89 mHit;
90 mHit=(AliMUONHit*)MUON->NextHit())
91 {
92 Int_t nch = mHit->fChamber; // chamber number
93 Float_t x = mHit->fX; // x-pos of hit
94 Float_t y = mHit->fY; // y-pos
95 Int_t ip = mHit->fParticle;
96 if (ip != kMuonPlus && ip != kMuonMinus) continue;
97 if (nch != 1) continue;
98 Float_t dmin=1.e8;
99 Int_t imin=-1;
100//
101// Loop over clusters
102 Int_t npos=0;
103
104 for (Int_t i=0; i<nrawcl; i++) {
105 Float_t dx = xc[i]-x;
106 Float_t dy = yc[i]-y;
107 Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
108 if (dr<dmin) {
109 imin=i;
110 dmin=dr;
111 }
112 if (TMath::Abs(dy)< 0.05 && TMath::Abs(dx)< 0.2 ) npos++;
113 }
114 mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(imin);
115 Int_t track=mRaw->fTracks[1];
116 Float_t xrec=mRaw->fX[1];
117 Float_t yrec=mRaw->fY[0];
118 Float_t x_res=xrec-x;
119 Float_t y_res=yrec-y;
120 hresx->Fill(x_res,1.);
121 hresy->Fill(y_res,1.);
122 hpos->Fill(Float_t(npos), 1.);
123 hresym->Fill(y_res*1.e4,1.);
124 if (npos == 0) {
125 printf("No cluster %d %f %f %d\n",
126 itrack, x, y, nhit);
127 h21->Fill(x,y,1.);
128 }
129
130
131
132 } // hits
133 } // tracks
134 } // end if MUON
135 } // event loop
136 TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
137 c1->Divide(2,2);
138 c1->cd(1);
139 hresx->SetFillColor(42);
140 hresx->SetXTitle("xrec-x");
141 hresx->Draw();
142
143 c1->cd(2);
144 hresy->SetFillColor(42);
145 hresy->SetXTitle("yrec-y");
146 hresy->Draw();
147
148 c1->cd(3);
149 hpos->SetFillColor(42);
150 hpos->SetXTitle("Possibilities");
151 hpos->Draw();
152
153 c1->cd(4);
154 hresym->SetFillColor(42);
155 hresym->SetXTitle("yrec-y");
156 hresym->Fit("gaus");
157 hresym->Draw();
158
159 TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
160 c2->Divide(2,2);
161
162 c2->cd(1);
163 h21->SetFillColor(42);
164 h21->SetXTitle("x");
165 h21->SetYTitle("y");
166 h21->Draw();
167
168 c2->cd(2);
169 hmult->SetFillColor(42);
170 hmult->SetXTitle("multiplicity");
171 hmult->Draw();
172
173 c2->cd(3);
174 hchi1->SetFillColor(42);
175 hchi1->SetXTitle("chi2");
176 hchi1->Draw();
177
178 c2->cd(4);
179 hchi2->SetFillColor(42);
180 hchi2->SetXTitle("chi2");
181 hchi2->Draw();
182
183}
184
185
186